# The Interpolating Element-Free Galerkin Method for 2D Transient Heat Conduction Problems

Authors: Na Zhao; Hongping Ren
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

An interpolating element-free Galerkin (IEFG) method is presented for transient heat conduction problems. The shape function in the moving least-squares (MLS) approximation does not satisfy the property of Kronecker delta function, so an interpolating moving least-squares (IMLS) method is discussed; then combining the shape function constructed by the IMLS method and Galerkin weak form of the 2D transient heat conduction problems, the interpolating element-free Galerkin (IEFG) method for transient heat conduction problems is presented, and the corresponding formulae are obtained. The main advantage of this approach over the conventional meshless method is that essential boundary conditions can be applied directly. Numerical results show that the IEFG method has high computational accuracy.


---

## Body

## 1. Introduction

In recent years, meshless methods have been successfully developed and applied to solve a variety of science and engineering problems [1–8]. The meshless method, which is based on nodes with a minimum of meshing or no meshing at all, can solve many engineering problems that are not suited to conventional computational methods and has shown some advantages.Many kinds of meshless methods have been developed, such as smoothed particle hydrodynamics (SPH) [9], diffuse element method (DEM) [10], element-free Galerkin (EFG) method [11], reproducing kernel particle method (RKPM) [12], finite point method (FPM) [13], meshless local Petrov-Galerkin (MLPG) method [14], point collocation method (PCM) [15], radial basis functions (RBF) method [16], meshless finite element method (MFEM) [17], complex variable meshless method (CVMM) [18, 19], boundary node method (BNM) [20], local boundary integral equation (LBIE) method [21], boundary radial point interpolation method (BRPIM) [22], and boundary element-free method (BEFM) [23–25].The element-free Galerkin (EFG) method is the most important meshless method, the shape function in the EFG method is formed with moving least-squares (MLS) approximation, a disadvantage of the MLS approximation is that the final algebraic equations system is sometimes ill-conditioned, and we cannot obtain a good solution, or even correctly obtain a numerical solution; then Cheng and Peng proposed an improved moving least-squares approximation by orthogonalizing the basis functions in the MLS approximation, and based on it Cheng and Peng put forward a boundary element-free method [23]. Since the shape function of the MLS approximation does not have the properties of Kronecker delta function, the meshless method based on it must use other methods, such as the penalty function method, Lagrange multiplier, to impose essential boundary conditions, which makes the weak form of the problem more complicated and the computational efficiency lower as a result, and Lancaster proposed interpolating moving least-squares (IMLS) method [26], which can obtain the shape function satisfying the property of Kronecker delta function, but, compared with the IMLS method, the shape function of the MLS approximation is much simpler, and thus a few papers on meshless methods based on the IMLS method were published. Ren et al. proposed an improved IMLS method, and based on it the interpolating element-free Galerkin (IEFG) method and the improved boundary element-free method are presented [27–30].The analysis of transient heat conduction problems is very important to engineering and science. However, analytical solution for this kind of problems is difficult to obtain except for a few simple cases; thus an alternative way is proposed, that is, the numerical solution. Singh et al. analyzed transient nonlinear heat transfer problems in solids by EFG method [31]. Chen and Cheng used complex variable reproducing kernel particle method to solve transient heat conduction problems [32]; its advantage is that 2D problem is solved with 1D basis function. Yang and Gao used radial integration BEM to solve transient heat conduction problems [33]; the features are that thermal material parameters can be functions of spatial coordinates. R. J. Cheng and Y. M. Cheng solved the inverse heat conduction problem [34]; the finite point method is used to obtain the solution of 1D inverse heat conduction problem with a source parameter. Chen and Liew, Baodong and Zheng developed meshless local Petrov-Galerkin approach based on the moving Kriging interpolation for transient heat conduction problems [35, 36]; the main feature is in use of the moving Kriging interpolation as the trail function and the Heaviside step function as the test function. Li et al. combined the complex variable reproducing kernel particle method and the finite element method to solve transient heat conduction problems [37]. Li et al. proposed the MLPG method in conjunction with the modified precise time step integration method for the analysis of transient heat conduction problems [38]. Zhang et al. developed improved element-free Galerkin method for three-dimensional transient heat conduction problems [39]; the effects of scaling parameter, number of nodes, and the time step length were considered in their numerical solutions.The present paper is motivated by the MLS method and the IMLS method; this is the fundamental principle of the EFG method and the IEFG method; then combining the shape function constructed by the IMLS method and Galerkin weak form of the 2D transient heat conduction problem, the interpolating element-free Galerkin (IEFG) method for 2D transient heat conduction problems is presented, and the corresponding formulae are obtained. For the purposes of demonstration, some selected numerical examples are solved using the IEFG method.


## 2. The Interpolating Moving Least-Squares (IMLS) Method

In the MLS approximation, it is assumed that a functionu
(
x
)
(
x
∈
D
) is to be approximated and that its values u
i
=
u
(
x
i
)
(
i
=
1,2
,
…
,
n
) are given.An approximating function ofu
(
x
) is
(1)
u
h
(
x
)
=
∑
i
=
1
m
p
i
(
x
)
a
i
(
x
)
=
p
T
(
x
)
a
(
x
)
,
where p
i
(
x
) are monomial basis functions and a
i
(
x
) are the coefficients of the basis functions (
i
=
1,2
,
…
,
m
).In general, the basis functions are as follows in 2D space.Linear basis:(2)
p
T
(
x
)
=
(
1
,
x
,
y
)
.Quadratic basis:(3)
p
T
(
x
)
=
(
1
,
x
,
y
,
x
2
,
x
y
,
y
2
)
.Letw
(
∥
x
-
x
I
∥
)
=
w
(
d
I
) be a weight function with compact support, and define a functional
(4)
J
=
∑
I
=
1
n
w
(
d
I
)
[
∑
i
=
1
m
p
i
(
x
I
)
·
a
i
(
x
)
-
u
(
x
I
)
]
2
,
and x
I
(
I
=
1,2
,
…
,
n
) are the nodes in the influence domain of point x and J is the sum of the squares of the residuals of all data points in the influence domain.WhenJ
=
min
⁡, we can obtain the coefficients a
i
(
x
):
(5)
a
(
x
)
=
A
-
1
(
x
)
B
(
x
)
u
,
where matrices A
(
x
) and B
(
x
) are
(6)
A
(
x
)
=
∑
I
=
1
n
w
(
d
I
)
p
(
x
I
)
p
T
(
x
I
)
=
P
T
W
(
x
)
P
,
B
(
x
)
=
(
w
(
d
1
)
p
(
x
1
)
,
w
(
d
2
)
p
(
x
3
)
,
…
,
w
(
d
n
)
p
(
x
n
)
)
B
(
x
)
=
P
T
W
(
x
)
,
p
(
x
)
=
(
p
1
(
x
)
,
p
2
(
x
)
,
…
,
p
m
(
x
)
)
T
.The expression of the approximation functionu
h
(
x
) is then
(7)
u
h
(
x
)
=
Φ
(
x
)
u
=
∑
I
=
1
n
Φ
I
(
x
)
u
I
,
where Φ
(
x
) is called the shape function and
(8)
Φ
(
x
)
=
(
Φ
1
(
x
)
,
Φ
2
(
x
)
,
…
,
Φ
n
(
x
)
)
=
p
T
(
x
)
A
-
1
(
x
)
B
(
x
)
.The normalized weight function is(9)
v
(
d
i
)
=
w
(
d
i
)
∑
I
=
1
n
w
(
d
I
)
(
i
=
1,2
,
…
,
n
)
.
Let
(10)
v
(
x
)
=
(
v
(
d
1
)
,
v
(
d
2
)
,
…
,
v
(
d
n
)
)
T
.
Assume
(11)
u
s
(
x
)
=
∑
I
=
1
n
u
(
x
I
)
v
(
d
I
)
=
v
T
(
x
)
u
,
which is a weighted average of the function values at nodes x
I in the influence domain of x.In the moving least-squares method, approximating functionu
h
(
x
) need not interpolate the data points; then by orthogonalizing the last m
-
1 basis functions to the first one, and taking a singular weight function in the points, we establish the interpolating moving least-squares (IMLS) method.We define the following inner product of functionf
(
x
) and g
(
x
):
(12)
(
f
,
g
)
x
=
∑
I
=
1
n
w
(
d
I
)
f
(
x
I
)
g
(
x
I
)
.Normalizep
1
(
x
)
≡
1 at point x; we have
(13)
b
x
(
1
)
=
p
1
∥
p
1
∥
x
=
1
[
∑
I
=
1
n
w
(
d
I
)
]
1
/
2
,
for i
=
2
,
3
,
…
,
m; generate the basis functions orthogonal to b
x
(
1
) as follows:
(14)
b
x
(
i
)
(
x
)
=
p
i
(
x
)
-
(
p
i
,
b
x
(
1
)
)
x
b
x
(
1
)
(
x
)
b
x
(
i
)
(
x
)
=
p
i
(
x
)
-
∑
I
=
1
n
p
i
(
x
I
)
w
(
d
I
)
∑
I
=
1
n
w
(
d
I
)
b
x
(
i
)
(
x
)
=
p
i
(
x
)
-
p
i
s
(
x
)
.In the MLS approximation, while we use the new basis functionsb
x
(
1
)
(
x
)
,
b
x
(
2
)
(
x
)
,
…
,
b
x
(
m
)
(
x
), the corresponding approximating function is
(15)
u
h
(
x
)
=
(
u
,
b
x
(
1
)
)
x
b
x
(
1
)
(
x
)
+
∑
i
=
2
m
a
i
-
1
(
x
)
b
x
(
i
)
(
x
)
u
h
(
x
)
=
v
T
(
x
)
u
+
b
T
(
x
)
A
x
-
1
(
x
)
B
x
(
x
)
u
=
Φ
(
x
)
u
.
Then
(16)
Φ
(
x
)
=
(
Φ
1
(
x
)
,
Φ
1
(
x
)
,
…
,
Φ
1
(
x
)
)
Φ
(
x
)
=
v
T
(
x
)
+
b
T
(
x
)
A
x
-
1
(
x
)
B
x
(
x
)
,
(
Φ
I
(
x
)
)
,
j
=
v
(
d
I
)
,
j
+
b
T
,
j
A
x
-
1
(
B
x
)
I
(
Φ
I
(
x
)
)
,
j
=
+
b
T
(
A
x
-
1
)
,
j
(
B
x
)
I
+
b
T
A
x
-
1
(
B
x
)
I
,
j
,
where
(17)
(
B
x
)
I
,
j
=
w
(
d
I
)
,
j
b
x
(
x
I
)
+
w
(
d
I
)
b
x
(
x
I
)
,
j
,
(
A
x
-
1
)
,
j
=
-
A
x
-
1
(
A
x
)
,
j
A
x
-
1
,
(
A
x
)
,
j
=
∑
I
=
1
n
w
(
d
I
)
,
j
b
x
(
x
I
)
b
x
T
(
x
I
)
(
A
x
)
,
j
=
+
∑
I
=
1
n
w
(
d
I
)
b
x
(
x
I
)
,
j
b
x
T
(
x
I
)
(
A
x
)
,
j
=
+
∑
I
=
1
n
w
(
d
I
)
b
x
(
x
I
)
b
x
T
(
x
I
)
,
j
.In the IMLS method, the singular weight functionw
(
d
) is selected as
(18)
w
(
d
)
{
1
d
2
(
1
-
1
d
)
2
,
d
≤
1
,
0
,
d
>
1
,
where d
=
∥
x
-
x
I
∥
/
r and r is the radius of the influence domain.In [4], we have proved that the approximating function u
h
(
x
) in (15) can interpolate at all points x
I, I
=
1,2
,
…
,
n; that is, u
h
(
x
I
)
=
u
(
x
I
).


## 3. 2D Transient Heat Conduction Analysis by Interpolating Element-Free Galerkin (IEFG) Method

### 3.1. Governing Equation and Its Galerkin Weak Form

In general, the governing equation for two-dimensional transient heat conduction in isotropic solid body with spatially varying conductivity occupying a regionΩ can be written as
(19)
c
ρ
∂
T
∂
t
=
∂
∂
x
(
k
1
∂
T
∂
x
)
+
∂
∂
y
(
k
2
∂
T
∂
y
)
+
Q
,
(
x
,
y
)
∈
Ω
,
where T
=
T
(
x
,
t
)
=
T
(
x
,
y
,
t
) represents temperature, t is time, ρ
(
x
) is the density of material, c
(
x
) is the specific heat capacity, Q
(
x
,
t
) is heat generation rate, and k
1 and k
2 are thermal conductivities in x- and y-directions, respectively.The initial condition is(20)
T
(
x
,
y
,
0
)
=
T
0
,
in
Ω
,
the Dirichlet boundary condition is
(21)
T
=
T
-
,
on
⁡
Γ
1
,
the Neumann boundary condition is
(22)
k
1
∂
T
∂
x
n
x
+
k
2
∂
T
∂
y
n
y
=
q
-
,
on
⁡
Γ
2
,
and the Robin boundary condition is
(23)
k
1
∂
T
∂
x
n
x
+
k
2
∂
T
∂
y
n
y
=
h
(
T
f
-
T
)
,
on
⁡
Γ
3
,
where (
n
x
,
n
y
) is the unit outward normal to the boundary Γ, T
- and q
- are the prescribed temperature and the given heat fluxes on the corresponding boundaries, h is the convection heat transfer coefficient, and T
f is the environmental temperature.Consider(24)
Γ
=
Γ
1
∪
Γ
2
∪
Γ
3
.The weak form of (19), (22), and (23) is
(25)
∫
Ω
∇
w
(
k
1
∂
T
∂
x
+
k
2
∂
T
∂
y
)
d
Ω
+
∫
Ω
w
(
c
ρ
∂
T
∂
t
-
Q
)
d
Ω
-
∫
Γ
2
w
·
q
-
d
Γ
-
∫
Γ
3
w
h
(
T
f
-
T
)
d
Γ
=
0
.
The functional Π
(
T
) can be written as
(26)
Π
(
T
)
=
1
2
∫
Ω
[
k
1
(
∂
T
∂
x
)
2
+
(
k
2
∂
T
∂
y
)
2
]
d
Ω
Π
(
T
)
=
+
∫
Ω
T
(
c
ρ
∂
T
∂
t
-
Q
)
d
Ω
Π
(
T
)
=
-
∫
Γ
2
T
·
q
-
d
Γ
-
∫
Γ
3
h
(
T
T
f
-
T
2
2
)
d
Γ
.
Let δ
Π
=
0; then
(27)
∫
Ω
δ
T
·
c
ρ
·
∂
T
∂
t
d
Ω
+
∫
Ω
δ
(
L
T
)
T
k
~
(
L
T
)
d
Ω
-
∫
Ω
δ
T
·
Q
d
Ω
-
∫
Γ
2
δ
T
·
q
-
d
Γ
-
∫
Γ
3
δ
T
·
h
(
T
f
-
T
)
d
Γ
=
0
,
where δ is the variational operator, and L is a differential operator,
(28)
L
(
·
)
=
[
∂
∂
x
∂
∂
y
]
(
·
)
,
k
~
=
[
k
1
0
0
k
2
]
.


### 3.2. Discretization of the Weak Form

We employedn
t nodes in the domain Ω, and the union of its compact support domain Ω
I, I
=
1,2
,
…
,
n
t, must cover Ω.From the approximating function (14), the temperature T
(
x
,
t
) at an arbitrary field point x
=
(
x
,
y
) in the domain Ω at a given time t can be expressed as
(29)
T
=
T
(
x
,
t
)
=
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
Φ
(
x
)
T
,
where n is the number of nodes in which the compact support domains cover the field point x.Moreover(30)
Φ
(
x
)
=
(
Φ
1
(
x
)
,
Φ
2
(
x
)
,
…
,
Φ
n
(
x
)
)
,
(31)
T
=
(
T
1
(
t
)
,
T
2
(
t
)
,
…
,
T
n
(
t
)
)
T
,
(32)
∂
T
∂
t
=
∂
∂
t
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
∑
I
=
1
n
Φ
I
(
x
)
∂
T
I
(
t
)
∂
t
=
Φ
(
x
)
T
˙
,
(33)
T
˙
=
(
∂
T
1
(
t
)
∂
t
,
∂
T
2
(
t
)
∂
t
,
…
,
∂
T
n
(
t
)
∂
t
)
T
,
(34)
L
T
=
L
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
∑
I
=
1
n
[
∂
∂
x
∂
∂
y
]
Φ
I
(
x
)
T
I
(
t
)
=
B
(
x
)
T
,
where
(35)
Β
(
x
)
=
(
Β
1
(
x
)
,
Β
2
(
x
)
,
…
,
Β
n
(
x
)
)
,
B
I
(
x
)
=
[
Φ
I
,
x
(
x
)
Φ
I
,
y
(
x
)
]
.
Substituting (29), (32), and (34) into (27) yields
(36)
∫
Ω
δ
[
Φ
(
x
)
T
]
c
ρ
[
Φ
(
x
)
T
˙
]
d
Ω
+
∫
Ω
δ
[
B
(
x
)
T
]
T
k
~
[
B
(
x
)
T
]
d
Ω
-
∫
Ω
δ
[
Φ
(
x
)
T
]
·
Q
d
Ω
-
∫
Γ
2
δ
[
Φ
(
x
)
T
]
·
q
-
d
Γ
-
∫
Γ
3
δ
[
Φ
(
x
)
T
]
·
h
T
f
d
Γ
+
∫
Γ
3
δ
[
Φ
(
x
)
T
]
h
[
Φ
(
x
)
T
]
d
Γ
=
0
.
The first term of (36) can be written as
(37)
∫
Ω
δ
[
Φ
(
x
)
T
]
c
ρ
[
Φ
(
x
)
T
˙
]
d
Ω
=
δ
T
T
[
∫
Ω
Φ
T
(
x
)
·
c
ρ
·
Φ
(
x
)
d
Ω
]
·
T
˙
=
δ
T
T
·
C
·
T
˙
,
(38)
C
=
∫
Ω
Φ
T
(
x
)
·
c
ρ
·
Φ
(
x
)
d
Ω
.
C
=
[
C
I
J
] is a matrix of n
t
×
n
t, where
(39)
C
I
J
=
∫
Ω
Φ
I
(
x
)
·
c
ρ
·
Φ
J
(
x
)
d
Ω
.The second term of (36) can be written as
(40)
∫
Ω
δ
[
B
(
x
)
T
]
T
k
~
[
B
(
x
)
T
]
d
Ω
=
δ
T
T
[
∫
Ω
B
T
(
x
)
·
k
~
·
B
(
x
)
d
Ω
]
·
T
=
δ
T
T
·
K
·
T
,
where
(41)
K
=
∫
Ω
B
T
(
x
)
·
k
~
·
B
(
x
)
d
Ω
.
K
=
[
K
I
J
] is a matrix of n
t
×
n
t, where
(42)
K
I
J
=
∫
Ω
[
k
1
∂
Φ
I
∂
x
∂
Φ
J
∂
x
+
k
2
∂
Φ
I
∂
y
∂
Φ
J
∂
y
]
d
Ω
.The third term of (36) can be written as
(43)
∫
Ω
δ
[
Φ
(
x
)
T
]
·
Q
d
Ω
=
δ
T
T
∫
Ω
Φ
T
·
Q
d
Ω
=
δ
T
T
·
F
(
1
)
,
(44)
F
(
1
)
=
(
f
1
(
1
)
(
t
)
,
f
2
(
1
)
(
t
)
,
…
,
f
n
(
1
)
(
t
)
)
T
,
(45)
f
I
(
1
)
(
t
)
=
∫
Ω
Φ
I
(
x
)
·
Q
d
Ω
.The fourth term of (36) can be written as
(46)
∫
Γ
2
δ
[
Φ
(
x
)
T
]
·
q
-
d
Γ
=
δ
T
T
∫
Γ
2
Φ
T
·
q
-
d
Γ
=
δ
T
T
·
F
(
2
)
,
(47)
F
(
2
)
=
(
f
1
(
2
)
(
t
)
,
f
2
(
2
)
(
t
)
,
…
,
f
n
(
2
)
(
t
)
)
T
,
(48)
f
I
(
2
)
(
t
)
=
∫
Γ
2
Φ
I
(
x
)
·
q
-
d
Γ
.The fifth term of (36) can be written as
(49)
∫
Γ
3
δ
[
Φ
(
x
)
T
]
·
h
T
f
d
Γ
=
δ
T
T
∫
Γ
3
Φ
T
·
h
T
f
d
Γ
=
δ
T
T
·
F
(
3
)
,
(50)
F
(
3
)
=
(
f
1
(
3
)
(
t
)
,
f
2
(
3
)
(
t
)
,
…
,
f
n
(
3
)
(
t
)
)
T
,
(51)
f
I
(
3
)
(
t
)
=
∫
Γ
3
Φ
I
(
x
)
·
h
T
f
d
Γ
.The sixth term of (36) can be written as
(52)
∫
Γ
3
δ
[
Φ
(
x
)
T
]
h
[
Φ
(
x
)
T
]
d
Γ
=
δ
T
T
[
∫
Γ
3
Φ
T
(
x
)
h
Φ
(
x
)
d
Γ
]
T
=
δ
T
T
·
H
·
T
.
H
=
[
H
I
J
] is a matrix of n
t
×
n
t, where
(53)
H
I
J
=
∫
Γ
3
Φ
I
(
x
)
·
h
·
Φ
J
(
x
)
d
Γ
.Substituting (37), (40), (43), (46), (49), and (52) into (36) yields
(54)
δ
T
T
·
(
C
T
˙
+
K
T
+
H
T
-
F
(
1
)
-
F
(
2
)
-
F
(
3
)
)
=
0
.
Because of the arbitrariness of δ
T
T, we have
(55)
C
T
˙
+
(
K
+
H
)
T
=
F
,
where
(56)
F
=
F
(
1
)
+
F
(
2
)
+
F
(
3
)
.Now, (27) has been discretized to be the ordinary differential equation (55), in which the time is the only variable.The traditional two-point difference method is selected for the time discretization,(57)
θ
(
∂
T
∂
t
)
t
+
Δ
t
+
(
1
-
θ
)
(
∂
T
∂
t
)
t
=
T
t
+
Δ
t
-
T
t
Δ
t
(
0
≤
θ
≤
1
)
.
Selecting θ
=
0, then (55) can be written as
(58)
C
T
n
+
1
-
T
n
Δ
t
+
(
K
+
H
)
T
n
+
1
+
T
n
2
=
F
n
+
1
+
F
n
2
;
that is,
(59)
(
C
+
Δ
t
(
K
+
H
)
2
)
T
n
+
1
=
(
C
-
Δ
t
(
K
+
H
)
2
)
T
n
+
Δ
t
(
F
n
+
1
+
F
n
)
2
,
where
(60)
T
n
+
1
=
T
(
(
n
+
1
)
Δ
t
)
,
T
n
=
T
(
n
Δ
t
)
,
F
n
+
1
=
F
(
(
n
+
1
)
Δ
t
)
,
F
n
=
F
(
n
Δ
t
)
.Substituting the boundary condition (21) into (59) directly, we can obtain the unknowns at nodes by solving (59).Above all, the IEFG method is presented for two-dimensional transient heat conduction problems.


## 3.1. Governing Equation and Its Galerkin Weak Form

In general, the governing equation for two-dimensional transient heat conduction in isotropic solid body with spatially varying conductivity occupying a regionΩ can be written as
(19)
c
ρ
∂
T
∂
t
=
∂
∂
x
(
k
1
∂
T
∂
x
)
+
∂
∂
y
(
k
2
∂
T
∂
y
)
+
Q
,
(
x
,
y
)
∈
Ω
,
where T
=
T
(
x
,
t
)
=
T
(
x
,
y
,
t
) represents temperature, t is time, ρ
(
x
) is the density of material, c
(
x
) is the specific heat capacity, Q
(
x
,
t
) is heat generation rate, and k
1 and k
2 are thermal conductivities in x- and y-directions, respectively.The initial condition is(20)
T
(
x
,
y
,
0
)
=
T
0
,
in
Ω
,
the Dirichlet boundary condition is
(21)
T
=
T
-
,
on
⁡
Γ
1
,
the Neumann boundary condition is
(22)
k
1
∂
T
∂
x
n
x
+
k
2
∂
T
∂
y
n
y
=
q
-
,
on
⁡
Γ
2
,
and the Robin boundary condition is
(23)
k
1
∂
T
∂
x
n
x
+
k
2
∂
T
∂
y
n
y
=
h
(
T
f
-
T
)
,
on
⁡
Γ
3
,
where (
n
x
,
n
y
) is the unit outward normal to the boundary Γ, T
- and q
- are the prescribed temperature and the given heat fluxes on the corresponding boundaries, h is the convection heat transfer coefficient, and T
f is the environmental temperature.Consider(24)
Γ
=
Γ
1
∪
Γ
2
∪
Γ
3
.The weak form of (19), (22), and (23) is
(25)
∫
Ω
∇
w
(
k
1
∂
T
∂
x
+
k
2
∂
T
∂
y
)
d
Ω
+
∫
Ω
w
(
c
ρ
∂
T
∂
t
-
Q
)
d
Ω
-
∫
Γ
2
w
·
q
-
d
Γ
-
∫
Γ
3
w
h
(
T
f
-
T
)
d
Γ
=
0
.
The functional Π
(
T
) can be written as
(26)
Π
(
T
)
=
1
2
∫
Ω
[
k
1
(
∂
T
∂
x
)
2
+
(
k
2
∂
T
∂
y
)
2
]
d
Ω
Π
(
T
)
=
+
∫
Ω
T
(
c
ρ
∂
T
∂
t
-
Q
)
d
Ω
Π
(
T
)
=
-
∫
Γ
2
T
·
q
-
d
Γ
-
∫
Γ
3
h
(
T
T
f
-
T
2
2
)
d
Γ
.
Let δ
Π
=
0; then
(27)
∫
Ω
δ
T
·
c
ρ
·
∂
T
∂
t
d
Ω
+
∫
Ω
δ
(
L
T
)
T
k
~
(
L
T
)
d
Ω
-
∫
Ω
δ
T
·
Q
d
Ω
-
∫
Γ
2
δ
T
·
q
-
d
Γ
-
∫
Γ
3
δ
T
·
h
(
T
f
-
T
)
d
Γ
=
0
,
where δ is the variational operator, and L is a differential operator,
(28)
L
(
·
)
=
[
∂
∂
x
∂
∂
y
]
(
·
)
,
k
~
=
[
k
1
0
0
k
2
]
.


## 3.2. Discretization of the Weak Form

We employedn
t nodes in the domain Ω, and the union of its compact support domain Ω
I, I
=
1,2
,
…
,
n
t, must cover Ω.From the approximating function (14), the temperature T
(
x
,
t
) at an arbitrary field point x
=
(
x
,
y
) in the domain Ω at a given time t can be expressed as
(29)
T
=
T
(
x
,
t
)
=
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
Φ
(
x
)
T
,
where n is the number of nodes in which the compact support domains cover the field point x.Moreover(30)
Φ
(
x
)
=
(
Φ
1
(
x
)
,
Φ
2
(
x
)
,
…
,
Φ
n
(
x
)
)
,
(31)
T
=
(
T
1
(
t
)
,
T
2
(
t
)
,
…
,
T
n
(
t
)
)
T
,
(32)
∂
T
∂
t
=
∂
∂
t
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
∑
I
=
1
n
Φ
I
(
x
)
∂
T
I
(
t
)
∂
t
=
Φ
(
x
)
T
˙
,
(33)
T
˙
=
(
∂
T
1
(
t
)
∂
t
,
∂
T
2
(
t
)
∂
t
,
…
,
∂
T
n
(
t
)
∂
t
)
T
,
(34)
L
T
=
L
∑
I
=
1
n
Φ
I
(
x
)
T
I
(
t
)
=
∑
I
=
1
n
[
∂
∂
x
∂
∂
y
]
Φ
I
(
x
)
T
I
(
t
)
=
B
(
x
)
T
,
where
(35)
Β
(
x
)
=
(
Β
1
(
x
)
,
Β
2
(
x
)
,
…
,
Β
n
(
x
)
)
,
B
I
(
x
)
=
[
Φ
I
,
x
(
x
)
Φ
I
,
y
(
x
)
]
.
Substituting (29), (32), and (34) into (27) yields
(36)
∫
Ω
δ
[
Φ
(
x
)
T
]
c
ρ
[
Φ
(
x
)
T
˙
]
d
Ω
+
∫
Ω
δ
[
B
(
x
)
T
]
T
k
~
[
B
(
x
)
T
]
d
Ω
-
∫
Ω
δ
[
Φ
(
x
)
T
]
·
Q
d
Ω
-
∫
Γ
2
δ
[
Φ
(
x
)
T
]
·
q
-
d
Γ
-
∫
Γ
3
δ
[
Φ
(
x
)
T
]
·
h
T
f
d
Γ
+
∫
Γ
3
δ
[
Φ
(
x
)
T
]
h
[
Φ
(
x
)
T
]
d
Γ
=
0
.
The first term of (36) can be written as
(37)
∫
Ω
δ
[
Φ
(
x
)
T
]
c
ρ
[
Φ
(
x
)
T
˙
]
d
Ω
=
δ
T
T
[
∫
Ω
Φ
T
(
x
)
·
c
ρ
·
Φ
(
x
)
d
Ω
]
·
T
˙
=
δ
T
T
·
C
·
T
˙
,
(38)
C
=
∫
Ω
Φ
T
(
x
)
·
c
ρ
·
Φ
(
x
)
d
Ω
.
C
=
[
C
I
J
] is a matrix of n
t
×
n
t, where
(39)
C
I
J
=
∫
Ω
Φ
I
(
x
)
·
c
ρ
·
Φ
J
(
x
)
d
Ω
.The second term of (36) can be written as
(40)
∫
Ω
δ
[
B
(
x
)
T
]
T
k
~
[
B
(
x
)
T
]
d
Ω
=
δ
T
T
[
∫
Ω
B
T
(
x
)
·
k
~
·
B
(
x
)
d
Ω
]
·
T
=
δ
T
T
·
K
·
T
,
where
(41)
K
=
∫
Ω
B
T
(
x
)
·
k
~
·
B
(
x
)
d
Ω
.
K
=
[
K
I
J
] is a matrix of n
t
×
n
t, where
(42)
K
I
J
=
∫
Ω
[
k
1
∂
Φ
I
∂
x
∂
Φ
J
∂
x
+
k
2
∂
Φ
I
∂
y
∂
Φ
J
∂
y
]
d
Ω
.The third term of (36) can be written as
(43)
∫
Ω
δ
[
Φ
(
x
)
T
]
·
Q
d
Ω
=
δ
T
T
∫
Ω
Φ
T
·
Q
d
Ω
=
δ
T
T
·
F
(
1
)
,
(44)
F
(
1
)
=
(
f
1
(
1
)
(
t
)
,
f
2
(
1
)
(
t
)
,
…
,
f
n
(
1
)
(
t
)
)
T
,
(45)
f
I
(
1
)
(
t
)
=
∫
Ω
Φ
I
(
x
)
·
Q
d
Ω
.The fourth term of (36) can be written as
(46)
∫
Γ
2
δ
[
Φ
(
x
)
T
]
·
q
-
d
Γ
=
δ
T
T
∫
Γ
2
Φ
T
·
q
-
d
Γ
=
δ
T
T
·
F
(
2
)
,
(47)
F
(
2
)
=
(
f
1
(
2
)
(
t
)
,
f
2
(
2
)
(
t
)
,
…
,
f
n
(
2
)
(
t
)
)
T
,
(48)
f
I
(
2
)
(
t
)
=
∫
Γ
2
Φ
I
(
x
)
·
q
-
d
Γ
.The fifth term of (36) can be written as
(49)
∫
Γ
3
δ
[
Φ
(
x
)
T
]
·
h
T
f
d
Γ
=
δ
T
T
∫
Γ
3
Φ
T
·
h
T
f
d
Γ
=
δ
T
T
·
F
(
3
)
,
(50)
F
(
3
)
=
(
f
1
(
3
)
(
t
)
,
f
2
(
3
)
(
t
)
,
…
,
f
n
(
3
)
(
t
)
)
T
,
(51)
f
I
(
3
)
(
t
)
=
∫
Γ
3
Φ
I
(
x
)
·
h
T
f
d
Γ
.The sixth term of (36) can be written as
(52)
∫
Γ
3
δ
[
Φ
(
x
)
T
]
h
[
Φ
(
x
)
T
]
d
Γ
=
δ
T
T
[
∫
Γ
3
Φ
T
(
x
)
h
Φ
(
x
)
d
Γ
]
T
=
δ
T
T
·
H
·
T
.
H
=
[
H
I
J
] is a matrix of n
t
×
n
t, where
(53)
H
I
J
=
∫
Γ
3
Φ
I
(
x
)
·
h
·
Φ
J
(
x
)
d
Γ
.Substituting (37), (40), (43), (46), (49), and (52) into (36) yields
(54)
δ
T
T
·
(
C
T
˙
+
K
T
+
H
T
-
F
(
1
)
-
F
(
2
)
-
F
(
3
)
)
=
0
.
Because of the arbitrariness of δ
T
T, we have
(55)
C
T
˙
+
(
K
+
H
)
T
=
F
,
where
(56)
F
=
F
(
1
)
+
F
(
2
)
+
F
(
3
)
.Now, (27) has been discretized to be the ordinary differential equation (55), in which the time is the only variable.The traditional two-point difference method is selected for the time discretization,(57)
θ
(
∂
T
∂
t
)
t
+
Δ
t
+
(
1
-
θ
)
(
∂
T
∂
t
)
t
=
T
t
+
Δ
t
-
T
t
Δ
t
(
0
≤
θ
≤
1
)
.
Selecting θ
=
0, then (55) can be written as
(58)
C
T
n
+
1
-
T
n
Δ
t
+
(
K
+
H
)
T
n
+
1
+
T
n
2
=
F
n
+
1
+
F
n
2
;
that is,
(59)
(
C
+
Δ
t
(
K
+
H
)
2
)
T
n
+
1
=
(
C
-
Δ
t
(
K
+
H
)
2
)
T
n
+
Δ
t
(
F
n
+
1
+
F
n
)
2
,
where
(60)
T
n
+
1
=
T
(
(
n
+
1
)
Δ
t
)
,
T
n
=
T
(
n
Δ
t
)
,
F
n
+
1
=
F
(
(
n
+
1
)
Δ
t
)
,
F
n
=
F
(
n
Δ
t
)
.Substituting the boundary condition (21) into (59) directly, we can obtain the unknowns at nodes by solving (59).Above all, the IEFG method is presented for two-dimensional transient heat conduction problems.


## 4. Numerical Examples

We select four numerical examples to demonstrate the applicability of the IEFG method in transient heat conduction problem.For the purpose of convergence studies, the root-mean-square (RMS) error is defined as(61)
r
=
1
n
∑
i
=
1
n
(
u
i
num
-
u
i
exact
)
2
,
where n is the number of sample points.


### 4.1. The Transient Heat Conduction without Heat Generation

The first example considered is the 2D transient heat conduction in a square domainΩ, the thermal conductivities k
1
=
k
2
=
1
0
3 W/m·°C, the specific heat capacity c
=
1
0
3 J/kg·°C, and the density ρ
=
1
0
3
kg
/
m
3. Its governing equation is
(62)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
=
0
,
x
∈
[
0
,
π
]
,
y
∈
[
0
,
π
]
.The boundary condition is(63)
T
(
x
,
y
,
t
)
|
(
x
,
y
)
∈
Γ
=
0
.The initial temperature is(64)
T
(
x
,
y
,
0
)
=
10
sin
x
sin
y
.The analytical solution of this problem is(65)
T
(
x
,
y
,
t
)
=
10
e
-
2
t
sin
x
sin
y
.As shown in Figure1, regular distribution of 13
×
13 nodes is arranged in the domain Ω and time step is Δ
t
=
0.02
s in the computing process. The numerical results and the analytical solution of temperature at y
=
π
/
2 when t
=
0.5
s, t
=
1
s, t
=
1.5
s, and t
=
2
s are shown in Figure 2. In addition, the values of root-mean-square errors at different times when the time step is Δ
t
=
0.002
s are given in Figure 3, which shows that numerical solution converges as the time increases.Figure 1
Node arrangement.Figure 2
Temperature distribution aty
=
π
/
2.Figure 3
Variation of root-mean-square error with time.


### 4.2. The Transient Heat Conduction without Heat Generation

Consider(66)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
=
0
,
x
∈
[
0,1
]
,
y
∈
[
0,1
]
.The boundary conditions are(67)
T
(
0
,
y
,
t
)
=
e
-
2
t
sin
y
,
T
(
1
,
y
,
t
)
=
e
-
2
t
sin
(
1
+
y
)
,
T
(
x
,
0
,
t
)
=
e
-
2
t
sin
x
,
T
(
x
,
1
,
t
)
=
e
-
2
t
sin
(
1
+
x
)
.The initial condition is(68)
T
(
x
,
y
,
0
)
=
sin
(
x
+
y
)
.The analytical solution of this problem is(69)
T
(
x
,
y
,
t
)
=
sin
(
x
+
y
)
e
-
2
t
.The distribution of nodes in this example is the same as that in Figure1. The time step is chosen as Δ
t
=
0.001
s. The numerical results and the analytical solution of temperature at x
=
0.25, x
=
0.5, and x
=
0.75 when t
=
0.01
s are plotted in Figure 4; we can conclude that the numerical results are in good agreement with the analytical solutions.Figure 4
Temperature distribution att
=
0.01
s.


### 4.3. The Transient Heat Conduction with Heat Generation

The third example is the transient heat conduction equation with heat generation; consider(70)
∂
T
∂
t
=
∂
2
T
∂
x
2
+
∂
2
T
∂
y
2
-
2
,
x
∈
[
0
,
π
]
,
y
∈
[
0
,
π
]
.The boundary conditions are(71)
T
(
0
,
y
,
t
)
=
0
,
T
(
π
,
y
,
t
)
=
π
2
,
T
(
x
,
0
,
t
)
=
T
(
x
,
π
,
t
)
=
x
2
.The initial temperature is(72)
T
(
x
,
y
,
0
)
=
x
2
+
sin
x
sin
y
.The analytical solution of this problem is(73)
T
(
x
,
y
,
t
)
=
e
-
2
t
sin
x
sin
y
+
x
2
.As Figure5 shows, 21
×
21 nodes are distributed in Ω. The time step is Δ
t
=
0.01
s. The numerical results and the analytical solution at y
=
π
/
2 when t
=
0.1
s, t
=
0.3
s, and t
=
0.5
s are shown in Figure 6; it is obvious that the numerical results are in excellent agreement with the analytical solutions.Figure 5
Node distribution.Figure 6
Temperature distribution aty
=
π
/
2.


### 4.4. The Transient Heat Conduction with Heat Generation and Lateral Heat Loss

Consider(74)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
-
(
1
+
t
2
)
T
=
Q
(
x
,
y
)
,
h
h
h
h
h
h
h
h
h
h
h
h
.
h
i
x
∈
[
0,1
]
,
y
∈
[
0,1
]
,
where
(75)
Q
(
x
,
y
)
=
(
2
π
2
-
t
2
-
2
)
e
-
t
sin
(
π
x
)
cos
⁡
(
π
y
)
.The boundary conditions are:(76)
T
(
0
,
y
,
t
)
=
0
,
T
(
1
,
y
,
t
)
=
0
,
T
(
x
,
0
,
t
)
=
e
-
t
sin
(
π
x
)
,
T
(
x
,
1
,
t
)
=
-
e
-
t
sin
(
π
x
)
.The initial temperature is(77)
T
(
x
,
y
,
0
)
=
sin
(
π
x
)
cos
⁡
(
π
y
)
.The analytical solution of this problem is(78)
T
(
x
,
y
,
t
)
=
e
-
t
sin
(
π
x
)
cos
⁡
(
π
y
)
.As Figure7 shows, 11
×
11 nodes are distributed in the rectangular domain Ω. The time step is Δ
t
=
0.001
s. The numerical results and the analytical solution at x
=
0.5 and y
=
0.8 when t
=
0.1
s, t
=
0.3
s, and t
=
0.5
s are shown in Figures 8 and 9, respectively; it can be found that IEFG method works well for transient heat conduction problems with heat generation and even with lateral heat loss.Figure 7
Node distribution.Figure 8
Temperature distribution atx
=
0.5.Figure 9
Temperature distribution aty
=
0.8.


## 4.1. The Transient Heat Conduction without Heat Generation

The first example considered is the 2D transient heat conduction in a square domainΩ, the thermal conductivities k
1
=
k
2
=
1
0
3 W/m·°C, the specific heat capacity c
=
1
0
3 J/kg·°C, and the density ρ
=
1
0
3
kg
/
m
3. Its governing equation is
(62)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
=
0
,
x
∈
[
0
,
π
]
,
y
∈
[
0
,
π
]
.The boundary condition is(63)
T
(
x
,
y
,
t
)
|
(
x
,
y
)
∈
Γ
=
0
.The initial temperature is(64)
T
(
x
,
y
,
0
)
=
10
sin
x
sin
y
.The analytical solution of this problem is(65)
T
(
x
,
y
,
t
)
=
10
e
-
2
t
sin
x
sin
y
.As shown in Figure1, regular distribution of 13
×
13 nodes is arranged in the domain Ω and time step is Δ
t
=
0.02
s in the computing process. The numerical results and the analytical solution of temperature at y
=
π
/
2 when t
=
0.5
s, t
=
1
s, t
=
1.5
s, and t
=
2
s are shown in Figure 2. In addition, the values of root-mean-square errors at different times when the time step is Δ
t
=
0.002
s are given in Figure 3, which shows that numerical solution converges as the time increases.Figure 1
Node arrangement.Figure 2
Temperature distribution aty
=
π
/
2.Figure 3
Variation of root-mean-square error with time.


## 4.2. The Transient Heat Conduction without Heat Generation

Consider(66)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
=
0
,
x
∈
[
0,1
]
,
y
∈
[
0,1
]
.The boundary conditions are(67)
T
(
0
,
y
,
t
)
=
e
-
2
t
sin
y
,
T
(
1
,
y
,
t
)
=
e
-
2
t
sin
(
1
+
y
)
,
T
(
x
,
0
,
t
)
=
e
-
2
t
sin
x
,
T
(
x
,
1
,
t
)
=
e
-
2
t
sin
(
1
+
x
)
.The initial condition is(68)
T
(
x
,
y
,
0
)
=
sin
(
x
+
y
)
.The analytical solution of this problem is(69)
T
(
x
,
y
,
t
)
=
sin
(
x
+
y
)
e
-
2
t
.The distribution of nodes in this example is the same as that in Figure1. The time step is chosen as Δ
t
=
0.001
s. The numerical results and the analytical solution of temperature at x
=
0.25, x
=
0.5, and x
=
0.75 when t
=
0.01
s are plotted in Figure 4; we can conclude that the numerical results are in good agreement with the analytical solutions.Figure 4
Temperature distribution att
=
0.01
s.


## 4.3. The Transient Heat Conduction with Heat Generation

The third example is the transient heat conduction equation with heat generation; consider(70)
∂
T
∂
t
=
∂
2
T
∂
x
2
+
∂
2
T
∂
y
2
-
2
,
x
∈
[
0
,
π
]
,
y
∈
[
0
,
π
]
.The boundary conditions are(71)
T
(
0
,
y
,
t
)
=
0
,
T
(
π
,
y
,
t
)
=
π
2
,
T
(
x
,
0
,
t
)
=
T
(
x
,
π
,
t
)
=
x
2
.The initial temperature is(72)
T
(
x
,
y
,
0
)
=
x
2
+
sin
x
sin
y
.The analytical solution of this problem is(73)
T
(
x
,
y
,
t
)
=
e
-
2
t
sin
x
sin
y
+
x
2
.As Figure5 shows, 21
×
21 nodes are distributed in Ω. The time step is Δ
t
=
0.01
s. The numerical results and the analytical solution at y
=
π
/
2 when t
=
0.1
s, t
=
0.3
s, and t
=
0.5
s are shown in Figure 6; it is obvious that the numerical results are in excellent agreement with the analytical solutions.Figure 5
Node distribution.Figure 6
Temperature distribution aty
=
π
/
2.


## 4.4. The Transient Heat Conduction with Heat Generation and Lateral Heat Loss

Consider(74)
∂
T
∂
t
-
∂
2
T
∂
x
2
-
∂
2
T
∂
y
2
-
(
1
+
t
2
)
T
=
Q
(
x
,
y
)
,
h
h
h
h
h
h
h
h
h
h
h
h
.
h
i
x
∈
[
0,1
]
,
y
∈
[
0,1
]
,
where
(75)
Q
(
x
,
y
)
=
(
2
π
2
-
t
2
-
2
)
e
-
t
sin
(
π
x
)
cos
⁡
(
π
y
)
.The boundary conditions are:(76)
T
(
0
,
y
,
t
)
=
0
,
T
(
1
,
y
,
t
)
=
0
,
T
(
x
,
0
,
t
)
=
e
-
t
sin
(
π
x
)
,
T
(
x
,
1
,
t
)
=
-
e
-
t
sin
(
π
x
)
.The initial temperature is(77)
T
(
x
,
y
,
0
)
=
sin
(
π
x
)
cos
⁡
(
π
y
)
.The analytical solution of this problem is(78)
T
(
x
,
y
,
t
)
=
e
-
t
sin
(
π
x
)
cos
⁡
(
π
y
)
.As Figure7 shows, 11
×
11 nodes are distributed in the rectangular domain Ω. The time step is Δ
t
=
0.001
s. The numerical results and the analytical solution at x
=
0.5 and y
=
0.8 when t
=
0.1
s, t
=
0.3
s, and t
=
0.5
s are shown in Figures 8 and 9, respectively; it can be found that IEFG method works well for transient heat conduction problems with heat generation and even with lateral heat loss.Figure 7
Node distribution.Figure 8
Temperature distribution atx
=
0.5.Figure 9
Temperature distribution aty
=
0.8.


## 5. Conclusions

The present study is concentrated on the interpolating element-free Galerkin (IEFG) method for 2D transient heat conduction problems; compared with the conventional EFG method, the essential boundary conditions are applied naturally and directly in the IEFG method, and thus the IEFG method gives a greater computational efficiency. Numerical results show that the IEFG method has high computational accuracy.



---

*Source: 712834-2014-02-26.xml*