- Email: [email protected]
Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements
Accepted Manuscript
Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements H.H. Zhang , S.Y. Han , L.F. Fan PII: DOI: Reference:
S0307-904X(17)30192-0 10.1016/j.apm.2017.03.043 APM 11684
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
23 November 2016 12 March 2017 21 March 2017
Please cite this article as: H.H. Zhang , S.Y. Han , L.F. Fan , Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.03.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
The numerical manifold method is further developed to solve 2D transient heat conduction problems. The Wachspress polygonal element is extended to thermal analysis.
Regular hexagonal mathematical elements are adopted in numerical examples.
Good agreement between our results and existing solutions is achieved.
The advantages of the NMM in discretization and accuracy are demonstrated.
AC
CE
PT
ED
M
AN US
CR IP T
1
ACCEPTED MANUSCRIPT
Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements H.H. Zhang a,*, S.Y. Han a, L.F. Fan b,* a
School of Civil Engineering and Architecture, Nanchang Hangkong University, Nanchang, Jiangxi 330063, China College of Architecture and Civil Engineering, Beijing University of Technology,
CR IP T
b
Beijing 100084, China
* Corresponding author: [email protected] (H.H. Zhang); [email protected] (L.F. Fan)
AN US
ABSTRACT
Due to the use of dual cover systems, i.e., the mathematical cover and the physical cover, the numerical manifold method (NMM) is able to solve physical problems with boundary-inconsistent meshes. Meanwhile, n-gons (n>4) are very impressive attributing to their greater flexibility in discretization, less sensitivity to volumetric
M
and shear locking, and better suitability for complex microstructures simulation. In
ED
the present paper, the NMM, combined with Wachspress-type hexagonal elements, is developed to solve 2D transient heat conduction problems. Based on the governing
PT
equations, the NMM temperature approximation and the modified variational principle, the NMM discrete formulations are deduced. The solution strategy to
CE
time-dependent global equations and the spatial integration scheme are presented. The advantages of the proposed approach in both discretization and accuracy are
AC
demonstrated through several typical examples with increasing complexity. The extension of polygonal elements in unsteady thermal analysis within the NMM is realized.
Keywords: 2D transient heat conduction; Numerical manifold method; Polygonal elements; Temperature; Regular hexagonal elements 1. Introduction Most heat transport phenomena occurred in practical engineering are transient and the solution to such kind of problems is a classical and also very important topic. 2
ACCEPTED MANUSCRIPT
Analytical and numerical approaches are two major choices for the calculation of thermal fields. Commonly, analytical methods are limited to linear cases with simple geometry, regular boundary or initial conditions [1-3]. As alternatives, numerical tools are widely adopted for more general problems. Among the frequently used computational methods, the finite difference method (FDM), the finite element method (FEM), the finite volume method (FVM), the
CR IP T
boundary element method (BEM) and the meshless methods are representatives. Dai and Nassar [4] developed a Crank-Nicholson-type FDM to solve unstable heat transport equation at the microscale. Wang [5] combined the maximum principle for differential equations with the FDM to obtain the upper and lower solutions to
AN US
transient heat conduction problems. Hien and Kleiber [6] analyzed unsteady linear heat transfer problems with the stochastic FEM. Chen et al. [7] simulated the time-dependent heat transport in the formation of electric arc sprayed coating by the FEM. Using the lattice Boltzmann method and the FVM, Mishra and Roy [8] studied
M
the transient conduction and radiation heat transfer in both 1D and 2D geometries. Guillot and McCool [9] investigated the effect of numerical boundary conditions on
ED
local error and convergence of the FVM for unsteady heat conduction. Sutradhar et al. [10] applied the Laplace transform and the BEM to solve unstable heat transport
PT
problems in 3D homogeneous materials and functionally graded materials (FGMs). Abreu et al. [11] adopted the BEM to tackle 2D transient heat conduction problems in
CE
homogeneous and FGMs together with the convolution quadrature method. Li et al. [12] considered the unsteady heat transfer in homogeneous materials through the
AC
meshless Petrov-Galerkin method and the modified precise time integration scheme. Sophy et al. [13] proposed a space-time diffuse approximation meshless method to study time-dependent heat conduction problems with static and moving heat sources. Recent years, the numerical manifold method (NMM) [14, 15] gradually comes into our sights because of its unique dual cover systems (i.e., the mathematical cover and the physical cover), which lead to a unified framework for both continuous and discontinuous analysis. The dominant characteristics of the NMM lie in: (1) the mathematical cover may be inconsistent with both exterior and interior domain 3
ACCEPTED MANUSCRIPT
boundaries, which eases the meshing task to some extent; (2) local physical properties (e.g., the displacement or temperature jump across crack surfaces, the strain or heat flux discontinuity across material interfaces) can be manifested in essence or represented by incorporating associated local functions into the approximation. Since its emergence, the NMM has been applied and extended to solve a variety of
propagation [29-31] and fluid flow [32-34].
CR IP T
problems in different areas, e.g., crack analysis [16-25], rock stability [26-28], wave
At the same time, compared with the mostly used and well developed triangular and quadrilateral elements, n-sided (n>4) polygons are also very attractive in many aspects [35-37], e.g., (1) they can provide more flexibility in domain discretization
AN US
and mesh refinement; (2) they are less sensitive to shear and volumetric locking; (3) they are very suitable in complex microstructures modeling (e.g., in the simulation of random heterogeneous medias containing polygonal inclusions and voids); (4) they are more tolerant to large localized deformations. In the past decades, the polygonal
M
elements have been embedded in different numerical methods to solve various problems, e.g., in the FEM for linear elasticity [38] and finite elasticity [35], in the
ED
smoothed FEM for solid mechanics [39], in the extended FEM [40] and NMM [19] for fracture analysis, in the scaled boundary FEM for ealstodynamics [41] and in the
PT
discontinuous Galerkin method involving heat flow [42]. In this paper, to further enlarge the thermal element library and reveal the
CE
properties of the polygonal elements in the NMM background, the Wachspress polygonal elements [43], are further extended to handle 2D transient heat conduction
AC
problems within the platform of the NMM. To this end, the rest of the paper is organized as follows. In Section 2, the governing equations are provided. Section 3 addresses the details about the polygonal NMM for the concerned problems. To verify the proposed method, several numerical examples are tested in Section 4. Finally in Section 5, the concluding remarks are made. 2. Statement of the problem Consider transient heat conduction in a homogeneous isotropic planar domain 4
ACCEPTED MANUSCRIPT
bounded by the contour T
q shown in Fig. 1. The governing partial
differential equation for this kind of problem is [12]
T x, y, t 2T x, y, t 2T x, y, t c k k Q t x 2 y 2
(1)
where ρ, c and k are, respectively, the density, the specific heat at constant pressure
CR IP T
and the thermal conductivity of the material. denotes partial derivative. T, t and Q stand for, respectively, the temperature, the time and the heat source. (x, y) are the Cartesian coordinates plotted in Fig. 1. q
n
AN US
q
y
x
o
M
T
T
ED
Fig. 1. Transient heat conduction in a 2D domain
The associated boundary conditions are
T T nx k ny q on q x y
(3)
CE
k
(2)
PT
T T on T
AC
where T and q are, respectively, the prescribed temperature on Dirichlet boundary
T and the enforced heat flux on Neumann boundary q . (nx , ny ) n is the
outward unit normal to the domain as illustrated in Fig. 1. As a time-dependent problem, the initial condition should also be provided, i.e.,
T x, y,0 T0 ( x, y)
(4)
3. Unsteady heat conduction simulation by the polygonal NMM 3.1 A sketch of the NMM 5
ACCEPTED MANUSCRIPT
The NMM was originally brought out by Shi [14] to study rock engineering problems with satisfactory efficiency and accuracy. The predominance of the NMM, versus other frequently used computational methods (e.g., the FEM and the BEM) results from its dual cover system, i.e., the mathematical cover (MC) and the physical cover (PC). To discretize a physical domain by the NMM, the MC is firstly generated by a series of mathematical patches (MPs) composed of user-defined mathematical
CR IP T
elements. Different from the finite element (FE) mesh, the MC may be inconsistent with the boundary, as long as it is large enough to cover the problem domain. Then, the physical patches (PPs) are obtained by the intersection of MPs and the physical region. All the PPs form the PC. Finally, the manifold elements (MEs) are produced
AN US
through the common zone of PPs.
To further demonstrate the above process, the discretization of the domain
1
2 in Fig. 2a is presented. To cover the whole region, the MC made up of
two hexagonal MPs, that is, M1 and M2 in Fig. 2b, is chosen. The intersection of the
M
physical domain and the MPs (see Fig. 2c) gives 4 PPs in Fig. 2d, that is, P1 and P2 from M1 and , P3 and P4 from M2 and . Above 4 PPs finally generate 6 MEs,
ED
labeled by E1 (P1), E2 (P1, P3), E3 (P3), E4 (P2), E5 (P2, P4) and E6 (P4) in Fig. 2e, where
1
M1
M2
CE
2
PT
the bracketed contents represent the connected PPs.
(b)
AC
(a)
(c)
P1
P3
P2
P4
(d) 6
ACCEPTED MANUSCRIPT
P1
P3
E1
E2
E3
P2
P4
E4
E5
E6
(e) Fig. 2. Illustration of the NMM discretization: (a) the physical domain, (b) the MC, (c) intersection of the physical domain and MPs, (d) generation of the PPs and (e) formations of the MEs
CR IP T
Based on the above concepts, the NMM field interpolation function can be obtained. For the present problem, the temperature in a certain ME E is approximately written as n
T h ( x, y, t ) wi ( x, y )T i ( x, y, t )
(5)
AN US
i 1
where n is the amount of PPs shared by E. wi ( x, y) is the weight function on the ith PP, inherited from that defined on the MP MI containing the very PP, i.e., wI ( x, y) and will be further discussed in the coming subsection. T i ( x, y, t ) is the local
M
function defined on the ith PP and is commonly chosen as
(6)
ED
T i ( x, y, t ) P( x, y)ai ( x, y, t )
where a i is the vector of the generalized degrees of freedom (DOFs) associated with
PT
the ith PP. P( x, y) is the polynomial basis matrix, i.e., P( x, y) 1, x, y, x 2 , xy, y 2 , ...
CE
(7)
3.2 Weight functions for the polygonal elements
AC
According to the NMM fundamentals [14], the weight function wI ( x, y) on the
MP MI should meet the following conditions, that is,
wI ( x, y) C 0 , ( x, y) M I
(8)
wI ( x, y) 0, ( x, y) M I
(9)
where C 0 is the set of all continuous functions.
7
ACCEPTED MANUSCRIPT
Besides, to assure the conforming property, the weight functions on all the MPs containing the point ( x, y) should still satisfy
wJ ( x, y ) 1
(10)
J ( x , y )M J
Although, the shape of mathematical elements may theoretically be user-defined,
CR IP T
to simplify the generation of MC and take advantage of the well-constructed FEs, the noted triangular and quadrilateral elements (or shape functions) are extensively adopted as the NMM mathematical elements (or weight functions) (e.g., see [14, 17, 23]). In the present, the polygonal FEs are further borrowed to solve time-dependent heat transfer problems in the NMM. The birth of the polygonal FEs in the 1970s owes
AN US
to the pioneer, Eugene L. Wachspress and the original intension is to improve the discretization flexibility [43], besides which, more advantages of polygonal elements have been gradually revealed as mentioned in Section 1. In recent years, after the work of Wachspress, several other polygonal FE shape functions have been proposed,
M
including the mean value coordinates by Floater [44], the maximum entropy approximation by Sukumar [45], the natural neighbor shape functions by Sukumar
ED
and Tabarraei [37] and the metric coordinates by Malsch et al. [46]. The comparisons between different polygonal FE shape functions were comprehensively performed by
PT
Sukumar and Malsch [47].
Basically, the above mentioned polygonal shape functions all meet the conditions
CE
listed in Eqs. (8)-(10). In view that it is the earliest and also the simplest form on convex elements [47], the rational Wachspress FE shape function is adopted as the
AC
NMM weight function. For the N-sided convex polygonal element in Fig. 3, the corresponding weight function at point p ( x, y) associated with the vertex vI is
wI ( x, y )
I ( x, y )
(11)
N
J 1
J
( x, y )
where
8
ACCEPTED MANUSCRIPT
I x, y
A vI 1 , vI , vI 1 cot I cot I 2 A vI 1 , vI , p A vI , vI 1 , p pv
(12)
I
cot I cot I
vI p vI vI 1
(13)
vI p vI vI 1 vI vI 1 vI p
(14)
vI vI 1 vI p
CR IP T
where the last expression in Eq. (12) was contributed by Meyer et al. [48].
A p1 , p2 , p3 denotes the area bounded by three vertices p1 , p2 and p3 . I is the included angle between the segments vI p and vI vI 1 , and I is that between
product and module.
AN US
vI vI 1 and vI p . “ ”, “ ” and “| |” are, respectively, the vector dot product, cross
N
1 2
N-1
ED
M
p( x, y)
I I
vI
vI 1
PT
Fig. 3. Illustration of an N-sided polygonal element
What should be further emphasized is that for the case when p lies on certain
CE
element edge, i.e., on side vI vI 1 or vI vI 1 , Eq. (12) is invalid. Under this
AC
circumstances, wI ( x, y) can be calculated through interpolation owing to its linear variation along element edges. 3.3 NMM global equations Presently, the NMM discrete equations for the concerned problems are derived by virtue of the modified variational principle [49]. The equivalent functional in correspondence with Eqs. (1) and (3) is 2 2 T 1 T 1 T (T ) cT k + k QT d qTd q t 2 x 2 y
9
(15)
ACCEPTED MANUSCRIPT
Further, seeing that the MC may not conform to the domain boundaries and the generalized DOFs in Eq. (6) may possess no physical meaning, we utilize the penalty method, adopted widely in the NMM literatures (e.g., see [14, 17, 26, 50]), to approximately impose the Dirichlet boundary condition in Eq. (2). Accordingly, we obtain the modified functional as
1 (T T ) (T T )d 2 T
(16)
CR IP T
* (T ) (T )
where is the penalty.
Based on the variational principle, only the true solution T can satisfy the extreme condition of * , i.e., * (T ) 0 , through which and Eqs. (5), (6) as well, the
AN US
NMM global formulations for transient heat conduction are finally deduced as KT + CT F
(17)
where T and T are, respectively, the vector of DOFs and their time derivatives. K, C and F represent, respectively, the thermal conductivity matrix, the heat capacity
M
matrix and the equivalent thermal load vector. They are firstly computed on individual ME and then summed. For some ME E,
K E E BT kBd
FE
E
(18)
T wi P c wi P d
wi P
T
Qd
PT
E
T wi P wi P d
ED
CE
TE
TE
(19)
T T wi P Td wi P qd E q
(20)
CE
where the superscript T symbolizes the matrix transpose. E is the domain occupied by E. TE and qE denote, respectively, the essential and natural boundary
AC
associated with E. The matrix B can be further expressed as
B B1 B2 ... Bi ... Bn
(21)
( wi P), x Bi ( wi P), y
(22)
with
3.4 Computer implementations
10
ACCEPTED MANUSCRIPT
To solve the linear system of the ordinary differential equations in Eq. (17), the commonly used Euler backward difference scheme [49] is adopted. Accordingly, supposing that the DOFs vector Tt at the instant t is known and the time step size is t , we can obtain Tt t at t t by
C C K t t Tt t Ft t Tt t t
CR IP T
(23)
To acquire K t t , C and Ft t in Eq. (23), the associated spatial integration scheme should also be provided. In the original NMM [14], simplex integration scheme is proposed to analytically compute the domain integration of polynomial
AN US
integrands. In view of that the Wachspress-type weight functions for n-gons (n>4) are rational, this scheme is infeasible herein. In this sense, the numerical integration strategy is used. Recent years, many studies have been carried out with regard to the numerical integration on polygonal elements, e.g., see [51-54]. In the present paper, a
M
conventional technique adopted in [40, 55] is taken. In the computation, every ME is firstly triangulated, and then regular Gaussian quadrature rules are applied on each
ED
triangle. Finally, the results on all the sub-elements add up to those on the associated ME.
PT
4. Numerical examples
In this section, to verify the feasibility and accuracy of the proposed method, four
CE
transient heat conduction cases, under Dirichlet boundary condition or mixed Dirichlet and Neumann boundary conditions, are conducted. As for the shape of the
AC
MP, theoretically, any arbitrarily shaped mathematical element is available to form the MPs, however, in light of that the MC can be inconsistent with the physical boundary and the accuracy of regular elements [56], the equilateral triangles and squares are mostly used in the NMM, e.g., as in [14, 17, 18, 23, 28]. In this work, on the ground that among all the regular n-gons (n>4), only hexagons can individually cover the plane, MPs composed of regular hexagonal mathematical elements are adopted throughout the analysis. Furthermore, to clearly quantify the discretization, besides the amount of PPs and MEs, the element size (the side length of a mathematical element) is also used. In addition, to avoid the rank deficiency trouble in the use of 11
ACCEPTED MANUSCRIPT higher-order local function [57, 58], the polynomial basis in Eq. (7) is chosen to be constant. 4.1 A square plate under Dirichlet boundary condition In the first example, heat conduction in a square plate without heat source is considered. As shown in Fig. 4, the width of the plate is L. The temperature T is enforced along the entire boundary and the initial condition is set such that T ( x, y,0) 10sin( x)sin( y) .
CR IP T
In the analysis, the following parameters are taken: L , 1.0, c 1.0 k 1.0 and T 0 . Accordingly, the theoretical temperature field is [12]
T ( x, y, t ) 10sin( x)sin( y)exp(2t )
(24)
the Fourier’s law as qx 10cos( x)sin( y ) exp(2t ) q y 10sin( x) cos( y ) exp(2t ) y
ρ, c, k
T
T
L
x
L
(a)
(b)
(c)
(d)
AC
CE
o
PT
A
ED
B
T
(25)
M
T
AN US
and the components of the heat flux vector, that is, qx and qy, can be derived through
Fig. 4. Transient heat conduction in a square plate under Dirichlet boundary condition: (a) physical model, (b) discretization when h=0.1, (c) discretization when h=0.2 and (d) discretization when h=0.4 12
ACCEPTED MANUSCRIPT
Firstly, the sensitivity study on the effect of the penalty number is investigated. For this purpose, the variations of temperature along the line y =π/2 are concerned with ∆t = 0.02, under the discretization in Fig. 4b (element size h = 0.1, 959 PPs and 429 MEs) at different λ, that is, λ =10k, 102k, 103k，104k，106k，108k and 1010k. The simulated temperatures when x =iπ/12 (i= 1, 2…6) at two instants, i.e., t =0.5 and 1.0, are, respectively, shown in Table 1 and 2, together with the exact solutions from Eq.
CR IP T
(24). Accordingly, we can find that the error of the present results with relatively small λ (e.g., 10k and 102k) is greater than that of larger λ (e.g., 106k and more) at a given point; besides, for the cases when λ>104k, the difference among NMM solutions is ignorable. On this ground, λ will be fixed at 1010k in the following computations.
λ/k =10
1.017469
1.909451
λ/k =102
0.970114
1.866591
3
0.965342
1.862248
λ/k =104
0.964864
1.861813
λ/k =106
0.964812
1.861766
8
0.964811
λ/k =1010
0.964811
Exact
0.952144
x=5π/12
x=π/2
2.676400
3.270243
3.634246
3.757009
2.638406
3.236455
3.603165
3.726853
2.634538
3.233001
3.599979
3.723760
2.634150
3.232655
3.599660
3.723450
2.634108
3.232617
3.599625
3.723416
1.861765
2.634107
3.232616
3.599624
3.723415
1.861765
2.634107
3.232616
3.599624
3.723415
1.839401
2.601305
3.185934
3.553445
3.678794
PT
λ/k =10
x=π/3
M
x=π/6
λ/k =10
x=π/4
ED
x=π/12
AN US
Table 1 Calculated temperature at selected points under different penalty values when t=0.5
Table 2 Calculated temperature at selected points under different penalty values when t=1.0 x=π/6
x=π/4
x=π/3
x=5π/12
x=π/2
λ/k =10
0.388076
0.728180
1.020464
1.246647
1.385214
1.431934
λ/k =102
0.363929
0.700224
0.989741
1.214063
1.351606
1.397997
3
0.361531
0.697431
0.986660
1.210788
1.348223
1.394579
λ/k =104
0.361291
0.697152
0.986352
1.210460
1.347884
1.394237
λ/k =106
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
8
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
10
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
0.350274
0.676678
0.956967
1.172040
1.307240
1.353353
CE
x=π/12
AC
λ/k =10
λ/k =10
λ/k =10 Exact
Then, the convergence study is performed concerning the effect of element size. Another two MCs with h=0.2 and 0.4 (see the discretization in Fig. 4c with 262 PPs 13
ACCEPTED MANUSCRIPT
and 110 MEs when h=0.2, and Fig. 4d with 89 PPs and 33 MEs when h=0.4) are further used under ∆t =0.02. The calculated temperatures at two sample points in Fig. 4a, that is, point A: (π/4, π/4) and point B: (π/2, π/2), are respectively, summarized in Fig. 5a and 5b, together with the exact ones. From these figures, we can find that with the decrease of h, the NMM results gradually converge to the analytical solutions. 10.0
Exact NMM-h=0.4 NMM-h=0.2 NMM-h=0.1
8.0 7.0
Temperature
Temperature
4.0
Exact NMM-h=0.4 NMM-h=0.2 NMM-h=0.1
9.0
CR IP T
5.0
3.0
2.0
6.0 5.0 4.0 3.0
1.0
2.0 1.0
0.0 1.0
2.0
3.0
4.0
AN US
0.0
0.0
5.0
Time
0.0
1.0
2.0
3.0
4.0
5.0
Time
(a) (b) Fig. 5. Computed temperatures of sample points at different instants and element size: (a) point A and (b) point B
Next, the effect of time step size on the solution accuracy is investigated. To this
M
end, three time steps, i.e., ∆t =0.1, 0.05 and 0.02, are sequentially examined. The
ED
discretization in Fig. 4b is used. The computed temperatures of point A and B at different instants are plotted in Fig. 6, from which the decay of temperature with time
PT
can be detected. For comparison, the exact results from Eq. (24) are also illustrated therein. Obviously, the temperatures at all time steps match well with the analytical
CE
ones; what’s more, with the decrease of time step, our results are getting closer to the exact values, which also conforms to the convergence rule of the backward difference
AC
scheme.
10.0
5.0
Exact NMM-time step= 0.10 NMM-time step= 0.05 NMM-time step= 0.02
8.0 7.0
Temperature
Temperature
4.0
3.0
2.0
Exact NMM-time step=0.10 NMM-time step=0.05 NMM-time step=0.02
9.0
6.0 5.0 4.0 3.0
1.0
2.0 1.0
0.0
0.0 0.0
1.0
2.0
3.0
4.0
0.0
5.0
1.0
2.0
3.0
Time
Time
14
4.0
5.0
ACCEPTED MANUSCRIPT (a) (b) Fig. 6. Computed temperatures of sample points at different instants and time steps: (a) point A and (b) point B
Lastly, the accuracy of heat flux is further tested under the discretization in Fig. 4b and ∆t = 0.02. The flux components at point A are plotted in Fig. 7, together with the exact ones by Eq. (25). Upon these figures, we can conclude that the present flux
CR IP T
solutions also conform well to the analytical values. -5.0
-5.0
Exact NMM
Exact NMM
-3.0
-3.0
qx
qy
-4.0
-4.0
-2.0
-1.0
-1.0
AN US
-2.0
0.0
0.0 0.0
1.0
2.0
3.0
4.0
0.0
5.0
Time
1.0
2.0
3.0
4.0
5.0
Time
Fig. 7. Computed heat flux of point A at different instants: (a) qx and (b) qy
4.2 A square plate under mixed boundary conditions
M
In this example, we consider a square plate of width L with zero initial
ED
temperature in Fig. 8a. Just as in Example 1, the heat source is also ignored. The temperature T is imposed along the left boundary and the heat flux q is prescribed
PT
on the bottom boundary with the rest boundaries insulated. Note that this problem was previously solved by Li et al. with the meshless local Petrov-Galerkin method (MLPG)
CE
[12].
AC
y
T
o
ρ, c, k
B
q
A
L
x
L
15
ACCEPTED MANUSCRIPT (a)
(b)
Fig. 8. Transient heat conduction in a square plate under mixed boundary conditions: (a) physical model and (b) the discretization when h=2.8
In the modeling, the associated parameters are chosen as: L 100, 1.0,
c 1.0, k 1000, T 0 and q 1000 . Mathematical element of size h=2.8 is applied to form the MPs and the discretization finally generates 1158 PPs and 525 MEs as displayed in Fig. 8b. For comparison, the same time step as in [12], i.e.,
CR IP T
t 0.05 is used. Two sample points with their coordinates: (100, 0) and (100, 50)
(denoted, respectively, as point A and B in Fig. 8a), are concerned. The computed temperatures when t [0, 2.0] for both points are, respectively, depicted in Fig. 9a and 9b, together with those under regular nodes in [12] (note that the MLPG results
AN US
were approximately extracted from Figs. 9 and 10 of reference [12] ). It is easy to achieve that our solutions agree well with the existing ones. 16.0
50.0
14.0
45.0 40.0
12.0
MLPG Present
10.0
25.0 20.0
10.0 5.0 0.0 -5.0 -0.2
0.0
0.2
0.4
0.6
0.8
ED
15.0
1.0
1.2
1.4
1.6
1.8
Temperature
30.0
M
Temperature
35.0
MLPG Present
8.0 6.0 4.0 2.0 0.0
2.0
-2.0 -0.2
2.2
0.0
0.2
0.4
0.6
0.8
Time
1.0
1.2
1.4
1.6
1.8
2.0
2.2
PT
Time
(a)
(b)
CE
Fig. 9. Computed temperatures of sample points at different instants of: (a) point A and (b) point B
4.3 A concave domain under time-dependent boundary condition
AC
Herein, a concave domain with heat source subjected to time-dependent essential
boundary condition is considered. The dimensions and material properties are denoted in Fig. 10a. The temperature T T (t ) is applied on the left boundary and other boundaries are insulated. As for the initial condition, T ( x, y,0) 0 .
16
ACCEPTED MANUSCRIPT y
5 4
ρ, c, k
3
7
Q
12 11
6
2
10
1
9
H
CR IP T
T (t )
13
8
x
W2
W1
AN US
(a)
W1
M
(b)
ED
Fig. 10. Transient heat conduction in a concave domain under time-dependent boundary condition: (a) physical model and (b) the discretization when h=0.02
In the simulation, W1 0.3 , W2 0.4 , H 0.6 , 1.0 , c 1.0, k 1.0 , Q 10 ,
PT
T (t ) 60t , t 0.02 and h =0.02. The discretized domain is given in Fig. 10b, which
involves 1280 PPs and 570 MEs. Temperatures at 13 sample points (see Table 3 for
CE
their coordinates) in Fig. 10a are tested and the computed results at t =1.0 and 10.0 are summarized in Table 3. Concerning that the same problems were tackled by Yao et al.
AC
with the precise integration boundary element method (PIBEM) [59], we also report their solutions in this table. As we can see, the present results are in nice consistence with the reference values and the maximum relative error (RE) is less than 0.53% for t =1.0 and 0.07% for t =10.0.
17
ACCEPTED MANUSCRIPT
Table 3 Calculated temperature at sample points by NMM and PIBEM temperature
sample point t =1.0
t =10.0
NMM
PIBEM
RE (%)
NMM
PIBEM
RE (%)
1: (0.15, 0.1)
55.7389
55.6649
0.133
595.3105
595.2820
0.005
2: (0.15, 0.2)
55.0830
55.0101
0.132
594.5542
594.5438
0.002
3: (0.15, 0.3)
54.4254
54.3731
0.096
593.7956
593.8255
0.005
4: (0.15, 0.4)
54.0200
53.9746
0.084
593.3275
593.3760
0.008
5: (0.15, 0.5)
53.8201
53.7779
0.078
593.0964
593.1536
0.010
6: (0.50, 0.3)
42.4838
42.3573
0.299
580.0919
580.2965
0.035
7: (0.50, 0.4)
42.4340
42.2967
0.325
580.0375
580.2317
0.033
8: (0.50, 0.5)
42.3966
42.2550
0.335
579.9968
580.1868
0.033
9: (0.85, 0.1)
34.6041
34.4424
0.469
570.8062
571.1672
0.063
10: (0.85, 0.2)
35.0822
34.9159
0.476
571.3773
571.7179
0.060
11: (0.85, 0.3)
35.5935
35.4115
0.514
571.9888
572.2947
0.053
12: (0.85, 0.4)
35.9529
35.7641
0.528
572.4200
572.7055
0.050
13: (0.85, 0.5)
36.1552
35.9660
0.526
572.6634
572.9410
0.048
ED
AN US
CR IP T
(x, y)
M
number:
4.4 A square domain with multiple holes
PT
In the last example, to further demonstrate the advantages of the proposed method in discretization and accuracy as well, a square domain of width W with five holes in
CE
Fig. 11a is studied. The side length of the diamond (square) hole is a while the circular holes are of the same radius b. Material parameters and other dimensions are
AC
all presented therein. The boundary conditions are prescribed such that the temperature along the left boundary is T , the heat flux q1 (t ) and q2 (t ) continuously enter the plate from the top side and the central hole boundaries, respectively, and the remaining edges are all insulated.
18
ACCEPTED MANUSCRIPT q1 (t )
ρ, c, k
W1
B
b
b y
o
C
T
q2
q2
a
b
W2
x
A
q2
W2
b
D
W1
Q(t) W1
W2
W2
W1
W
(a)
CR IP T
q2
(b)
AN US
Fig. 11. Transient heat conduction in a square domain with five holes: (a) physical model and (b) the MC when h=0.15
In the NMM computations, the corresponding parameters are chosen as:
W 10.0, W1 W2 2.5 , a 1.5, b 1.0 , 1.0 , c 1.0,
k 1.0 , T 0 ,
q1 10(1 t ) and q2 5 . At the same time, the initial temperatures are set to be zero
M
and the heat source is Q(t ) 10sin(2t ) . Regular MC with h =0.15 as in Fig. 11b is
ED
used, regardless the size, location and boundaries of all interior holes and the exterior domain boundary as well. After some geometrical operations, the final NMM
AC
CE
PT
discretization is obtained in Fig. 12a, which includes 3457 PPs and 1567 MEs.
(a)
(b)
Fig. 12. Discretization by NMM and FEM: (a) NMM discretization and (b) FEM mesh
19
ACCEPTED MANUSCRIPT
To the authors’ knowledge, there are no reference solutions to this problem. Therefore, our results are compared with those from the widely used FE software ANSYS (version R14.5). The corresponding FE mesh is shown in Fig. 11b, with 2744 nodes and 2476 4-noded quadrilateral elements. From Figs. 11b, 12a and 12b, the unique feature of the NMM in meshing, as addressed in Section 1 and 3, can be easily inferred. In the comparative study, the time step in both methods is fixed at t 0.02 .
CR IP T
In addition, the Newmark scheme with ANSYS default parameters is used in the FEM analysis. The temperatures at 4 sample points in Fig. 11a, i.e., point A: (2.5, 0.0), point B: (0.0, 2.5), point C: (-2.5, 0.0) and point D: (0.0, -2.5), are, respectively, examined. The calculated results are plotted in Fig. 13a~d. As we can see, the
AN US
difference of temperatures by the NMM and the FEM is almost indistinguishable at all points.
45
18
FEM NMM
16
35
14 12 10 8 6
2 0
0.0
PT
-2 1.0
2.0
3.0
4.0
25 20 15 10
ED
4
Temperature
30
M
Temperature
FEM NMM
40
5 0 -5 5.0
0.0
1.0
2.0
Time
(a)
CE
16
4.0
5.0
(b) 16
FEM NMM
14
FEM NMM
14
12
12
10
10
Temperature
AC Temperature
3.0
Time
8 6 4
8 6 4
2
2
0
0
-2 0.0
1.0
2.0
3.0
4.0
-2
5.0
0.0
Time
1.0
2.0
3.0
Time
(c)
(d)
20
4.0
5.0
ACCEPTED MANUSCRIPT Fig. 13. The calculated temperatures of sample points at different instants by NMM and FEM: (a) point A, (b) point B, (c) point C and (d) point D
5. Concluding remarks In this work, the numerical manifold method was developed to tackle 2D unsteady heat transfer problems on hexagonal elements. Besides the description of governing equations, the introduction of the NMM, and the derivation of discrete formulations,
CR IP T
details about the weight functions of the Wachspress polygonal elements, the solution strategy to the time-dependent system of equations and the numerical integration schemes were also provided. Four numerical examples with increasing complexity are tested, under Dirichlet boundary condition or mixed Dirichlet and Neumann boundary conditions. Mathematical covers composed of regular hexagonal mathematical
AN US
elements were adopted throughout the simulation. The effect of some key parameters, i.e., the penalty value, the element size and the time step, on the results was investigated. The computed thermal fields agree well with the available reference solutions.
M
For the considered continuous heat transfer problems, compared with some other represented numerical tools, e.g., the FEM, the BEM and the meshless methods, the
ED
dominant advantage of the proposed approach lies in discretization, that is, the mathematical cover can be inconsistent with both external and internal boundaries
PT
(e.g., see example 4), which may reduce the meshing cost in part and exploit the accuracy of regular elements as well. Although only regular hexagons were used in
CE
this study, all convex n-gons (n>4) and their combinations are welcome from the essence of the Wachspress type polygonal NMM. The work on arbitrary polygonal
AC
elements is proceeding and the results are awaited. Acknowledgements The present work was supported by the National Natural Science Foundation of China (Grant Nos. 11462014 and 11572282), the Provincial Natural Science Foundation of Jiangxi, China (Grant No. 20151BAB202003), the Science and Technology Program of Educational Committee of Jiangxi Province of China (Grant Nos. GJJ14526 and GJJ150752). 21
ACCEPTED MANUSCRIPT Reference [1] M.N. Ozisik, Boundary values problems of heat conduction, Internatinal Textbook Company, Scranton, 1968. [2] S. Singh, P.K. Jain, Rizwan-uddin, Analytical solution to transient heat conduction in polar coordinates with multiple layers in radial direction, Int J Therm Sci, 47 (2008) 261-273. [3] H.J. Jiang, H.L. Dai, Analytical solutions for three-dimensional steady and transient heat conduction problems of a double-layer plate with a local heat source, Int J Heat Mass Tran, 89 (2015) 652-666. [4] W.Z. Dai, R. Nassar, A finite difference scheme for solving the heat transport equation at the
CR IP T
microscale, Numer Meth Part D E, 15 (1999) 697-708.
[5] C.C. Wang, Application of the maximum principle for differential equations in combination with the finite difference method to find transient approximate solutions of heat equations and error analysis, Numer Heat Tr B-Fund, 55 (2009) 56-72.
[6] T.D. Hien, M. Kleiber, Stochastic finite element modelling in linear transient heat transfer, Comput Method Appl M, 144 (1997) 111-124.
AN US
[7] Y.X. Chen, X.B. Liang, Y. Liu, J.Y. Bai, B.S. Xu, Finite element modeling of coating formation and transient heat transfer in the electric arc spray process, Int J Heat Mass Tran, 53 (2010) 2012-2021. [8] S.C. Mishra, H.K. Roy, Solving transient conduction and radiation heat transfer problems using the lattice Boltzmann method and the finite volume method, J Comput Phys, 223 (2007) 89-107. [9] M.J. Guillot, S.C. McCool, Effect of boundary condition approximation on convergence and accuracy of a finite volume discretization of the transient heat conduction equation, Int J Numer Method H, 25 (2015) 950-972.
M
[10] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method, Eng Anal
ED
Bound Elem, 26 (2002) 119-132.
[11] A.I. Abreu, A. Canelas, W.J. Mansur, A CQM-based BEM for transient heat conduction problems in homogeneous materials and FGMs, Appl. Math. Model., 37 (2013) 776-792.
PT
[12] Q.H. Li, S.S. Chen, G.X. Kou, Transient heat conduction analysis using the MLPG method and modified precise time step integration method, J Comput Phys, 230 (2011) 2736-2750. [13] T. Sophy, A. Da Silva, A. Kribeche, H. Sadat, An alternative space-time meshless method for
CE
solving transient heat transfer problems with high discontinuous moving sources, Numer Heat Tr B-Fund, 69 (2016) 377-388. [14] G.H. Shi, Manifold method of material analysis,
Transcations of the 9th Army Confernece on
AC
Applied Mathematics and ComputingMinneapolis, Minnesota, 1991, pp. 57-76. [15] G.W. Ma, X.M. An, L. He, The numerical manifold method: a review, Int J Comp Meth-Sing, 7 (2010) 1-32. [16] R.J. Tsay, Y.J. Chiou, W.L. Chuang, Crack growth prediction by manifold method, J Eng Mech-Asce, 125 (1999) 884-890. [17] G.W. Ma, X.M. An, H.H. Zhang, L.X. Li, Modeling complex crack problems using the numerical manifold method, International Journal of Fracture, 156 (2009) 21-35. [18] H.H. Zhang, L.X. Li, X.M. An, G.W. Ma, Numerical analysis of 2-D crack propagation problems using the numerical manifold method, Eng Anal Bound Elem, 34 (2010) 41-50. [19] H.H. Zhang, S.Q. Zhang, Extract of stress intensity factors on honeycomb elements by the numerical manifold method, Finite Elem Anal Des, 59 (2012) 55-65. 22
ACCEPTED MANUSCRIPT
[20] X.M. An, Z.Y. Zhao, H.H. Zhang, L. He, Modeling bimaterial interface cracks using the numerical manifold method, Eng Anal Bound Elem, 37 (2013) 464-474. [21] H. Zheng, D.D. Xu, New strategies for some issues of numerical manifold method in simulation of crack propagation, Int J Numer Meth Eng, 97 (2014) 986-1010. [22] H.H. Zhang, G.W. Ma, F. Ren, Implementation of the numerical manifold method for thermo-mechanical fracture of planar solids, Eng Anal Bound Elem, 44 (2014) 45-54. [23] H.H. Zhang, G.W. Ma, Fracture modeling of isotropic functionally graded materials by the numerical manifold method, Eng Anal Bound Elem, 38 (2014) 61-71. [24] J. He, Q.S. Liu, G.W. Ma, B. Zeng, An improved numerical manifold method incorporating hybrid
CR IP T
crack element for crack propagation simulation, International Journal of Fracture, 199 (2016) 21-38.
[25] H. Zheng, F. Liu, X.L. Du, Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method, Comput Method Appl M, 295 (2015) 150-171.
[26] G.X. Zhang, Y. Zhao, X.C. Peng, Simulation of toppling failure of rock slope by numerical manifold method, Int J Comp Meth-Sing, 7 (2010) 167-189.
[27] Y.J. Ning, X.M. An, G.W. Ma, Footwall slope stability analysis with the numerical manifold
AN US
method, Int J Rock Mech Min, 48 (2011) 964-975.
[28] L.N.Y. Wong, Z.J. Wu, Application of the numerical manifold method to model progressive failure in rock slopes, Eng Fract Mech, 119 (2014) 1-20.
[29] L.F. Fan, X.W. Yi, G.W. Ma, Numerical manifold method (Nmm) simulation of stress wave propagation through fractured rock mass, Int J Appl Mech, 5 (2013) 238-249.
[30] G.F. Zhao, X.B. Zhao, J.B. Zhu, Application of the numerical manifold method for stress wave propagation across rock masses, Int J Numer Anal Met, 38 (2014) 92-110.
M
[31] X.L. Qu, Y. Wang, G.Y. Fu, G.W. Ma, Efficiency and accuracy verification of the explicit numerical manifold method for dynamic problems, Rock Mech Rock Eng, 48 (2015) 1131-1142. [32] Q.H. Jiang, S.S. Deng, C.B. Zhou, W.B. Lu, Modeling unconfined seepage flow using
ED
three-dimensional numerical manifold method, J Hydrodyn, 22 (2010) 554-561. [33] H. Zheng, F. Liu, C.G. Li, Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method, Appl. Math. Model., 39 (2015) 794-808.
PT
[34] Y. Wang, M.S. Hu, Q.L. Zhou, J. Rutqvist, A new second-order numerical manifold method model with an efficient scheme for analyzing free surface flow with inner drains, Appl. Math. Model., 40
CE
(2016) 1427-1445.
[35] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, Polygonal finite elements for finite elasticity, Int J Numer Meth Eng, 101 (2015) 305-328.
AC
[36] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, Comput Method Appl M, 306 (2016) 216-251.
[37] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Int J Numer Meth Eng, 61 (2004) 2045-2066.
[38] A. Tabarraei, N. Sukumar, Application of polygonal finite elements in linear elasticity, Int J Comp Meth-Sing, 3 (2006) 503-520. [39] K.Y. Dai, G.R. Liu, T.T. Nguyen, An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics, Finite Elem Anal Des, 43 (2007) 847-860. [40] A. Tabarraei, N. Sukumar, Extended finite element method on polygonal and quadtree meshes, Comput Method Appl M, 197 (2008) 425-438. 23
ACCEPTED MANUSCRIPT
[41] Z.H. Zhang, Z.J. Yang, J.H. Li, An adaptive polygonal scaled boundary finite element method for elastodynamics, Int J Comp Meth-Sing, 13 (2016). [42] J. Jaskowiec, P. Plucinski, A. Stankiewicz, Discontinuous Galerkin method with arbitrary polygonal finite elements, Finite Elem Anal Des, 120 (2016) 1-17. [43] E.L. Wachspress, A rational finite element basis, Academic Press, New York, 1975. [44] M.S. Floater, Mean value coordinates, Comput Aided Geom D, 20 (2003) 19-27. [45] N. Sukumar, Construction of polygonal interpolants: a maximum entropy approach, Int J Numer Meth Eng, 61 (2004) 2159-2181. polygons, Journal of Graphics Tools, 10 (2005) 27-39.
CR IP T
[46] E.A. Malsch, J.J. Lin, G. Dasgupta, Smooth two-dimensional interpolations: a recipe for all [47] N. Sukumar, E.A. Malsch, Recent advances in the construction of polygonal finite element interpolants, Arch Comput Method E, 13 (2006) 129-163.
[48] M. Meyer, A. Barr, H. Lee, M. Desbrun, Generalized barycentric coordinates on irregular polygon, Journal of Graphics Tools, 7 (2002) 13-22.
[49] J.F. Wang, Y.M. Cheng, A new complex variable meshless method for transient heat conduction
AN US
problems, Chinese Phys B, 21 (2012).
[50] H.H. Zhang, G.W. Ma, L.F. Fan, Thermal shock analysis of 2D cracked solids using the numerical manifold method and precise time integration, Eng Anal Bound Elem, 75 (2017) 46-56. [51] G. Dasgupta, Integration within polygonal finite elements, J Aerospace Eng, 16 (2003) 9-18. [52] S. Natarajan, S. Bordas, D.R. Mahapatra, Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping, Int J Numer Meth Eng, 80 (2009) 103-134. [53] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on
M
irregular convex polygons and polyhedrons, Comput Mech, 47 (2011) 535-554. [54] E.B. Chin, J.B. Lasserre, N. Sukumar, Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra, Comput Mech, 56 (2015) 967-981.
ED
[55] A.R. Khoei, R. Yasbolaghi, S.O.R. Biabanaki, A polygonal finite element method for modeling crack propagation with minimum remeshing, International Journal of Fracture, 194 (2015) 123-148. [56] H.H. Zhang, Y.L. Chen, L.X. Li, X.M. An, G.W. Ma, Accuracy comparison of rectangular and
PT
triangular mathematical elements in the numerical manifold method, Analysis of Discontinuous Deformation: New Developments and Applications, (2010) 297-303.
CE
[57] R. Tian, G. Yagawa, H. Terasaka, Linear dependence problems of partition of unity-based generalized FEMs, Comput Method Appl M, 195 (2006) 4768-4782. [58] X.M. An, L.X. Li, G.W. Ma, H.H. Zhang, Prediction of rank deficiency in partition of unity-based
AC
methods with plane triangular or quadrilateral meshes, Comput Method Appl M, 200 (2011) 665-674. [59] W.A. Yao, B. Yu, X.W. Gao, Q. Gao, A precise integration boundary element method for solving transient heat conduction problems, Int J Heat Mass Tran, 78 (2014) 883-891.
24
Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements H.H. Zhang , S.Y. Han , L.F. Fan PII: DOI: Reference:
S0307-904X(17)30192-0 10.1016/j.apm.2017.03.043 APM 11684
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
23 November 2016 12 March 2017 21 March 2017
Please cite this article as: H.H. Zhang , S.Y. Han , L.F. Fan , Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.03.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
The numerical manifold method is further developed to solve 2D transient heat conduction problems. The Wachspress polygonal element is extended to thermal analysis.
Regular hexagonal mathematical elements are adopted in numerical examples.
Good agreement between our results and existing solutions is achieved.
The advantages of the NMM in discretization and accuracy are demonstrated.
AC
CE
PT
ED
M
AN US
CR IP T
1
ACCEPTED MANUSCRIPT
Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements H.H. Zhang a,*, S.Y. Han a, L.F. Fan b,* a
School of Civil Engineering and Architecture, Nanchang Hangkong University, Nanchang, Jiangxi 330063, China College of Architecture and Civil Engineering, Beijing University of Technology,
CR IP T
b
Beijing 100084, China
* Corresponding author: [email protected] (H.H. Zhang); [email protected] (L.F. Fan)
AN US
ABSTRACT
Due to the use of dual cover systems, i.e., the mathematical cover and the physical cover, the numerical manifold method (NMM) is able to solve physical problems with boundary-inconsistent meshes. Meanwhile, n-gons (n>4) are very impressive attributing to their greater flexibility in discretization, less sensitivity to volumetric
M
and shear locking, and better suitability for complex microstructures simulation. In
ED
the present paper, the NMM, combined with Wachspress-type hexagonal elements, is developed to solve 2D transient heat conduction problems. Based on the governing
PT
equations, the NMM temperature approximation and the modified variational principle, the NMM discrete formulations are deduced. The solution strategy to
CE
time-dependent global equations and the spatial integration scheme are presented. The advantages of the proposed approach in both discretization and accuracy are
AC
demonstrated through several typical examples with increasing complexity. The extension of polygonal elements in unsteady thermal analysis within the NMM is realized.
Keywords: 2D transient heat conduction; Numerical manifold method; Polygonal elements; Temperature; Regular hexagonal elements 1. Introduction Most heat transport phenomena occurred in practical engineering are transient and the solution to such kind of problems is a classical and also very important topic. 2
ACCEPTED MANUSCRIPT
Analytical and numerical approaches are two major choices for the calculation of thermal fields. Commonly, analytical methods are limited to linear cases with simple geometry, regular boundary or initial conditions [1-3]. As alternatives, numerical tools are widely adopted for more general problems. Among the frequently used computational methods, the finite difference method (FDM), the finite element method (FEM), the finite volume method (FVM), the
CR IP T
boundary element method (BEM) and the meshless methods are representatives. Dai and Nassar [4] developed a Crank-Nicholson-type FDM to solve unstable heat transport equation at the microscale. Wang [5] combined the maximum principle for differential equations with the FDM to obtain the upper and lower solutions to
AN US
transient heat conduction problems. Hien and Kleiber [6] analyzed unsteady linear heat transfer problems with the stochastic FEM. Chen et al. [7] simulated the time-dependent heat transport in the formation of electric arc sprayed coating by the FEM. Using the lattice Boltzmann method and the FVM, Mishra and Roy [8] studied
M
the transient conduction and radiation heat transfer in both 1D and 2D geometries. Guillot and McCool [9] investigated the effect of numerical boundary conditions on
ED
local error and convergence of the FVM for unsteady heat conduction. Sutradhar et al. [10] applied the Laplace transform and the BEM to solve unstable heat transport
PT
problems in 3D homogeneous materials and functionally graded materials (FGMs). Abreu et al. [11] adopted the BEM to tackle 2D transient heat conduction problems in
CE
homogeneous and FGMs together with the convolution quadrature method. Li et al. [12] considered the unsteady heat transfer in homogeneous materials through the
AC
meshless Petrov-Galerkin method and the modified precise time integration scheme. Sophy et al. [13] proposed a space-time diffuse approximation meshless method to study time-dependent heat conduction problems with static and moving heat sources. Recent years, the numerical manifold method (NMM) [14, 15] gradually comes into our sights because of its unique dual cover systems (i.e., the mathematical cover and the physical cover), which lead to a unified framework for both continuous and discontinuous analysis. The dominant characteristics of the NMM lie in: (1) the mathematical cover may be inconsistent with both exterior and interior domain 3
ACCEPTED MANUSCRIPT
boundaries, which eases the meshing task to some extent; (2) local physical properties (e.g., the displacement or temperature jump across crack surfaces, the strain or heat flux discontinuity across material interfaces) can be manifested in essence or represented by incorporating associated local functions into the approximation. Since its emergence, the NMM has been applied and extended to solve a variety of
propagation [29-31] and fluid flow [32-34].
CR IP T
problems in different areas, e.g., crack analysis [16-25], rock stability [26-28], wave
At the same time, compared with the mostly used and well developed triangular and quadrilateral elements, n-sided (n>4) polygons are also very attractive in many aspects [35-37], e.g., (1) they can provide more flexibility in domain discretization
AN US
and mesh refinement; (2) they are less sensitive to shear and volumetric locking; (3) they are very suitable in complex microstructures modeling (e.g., in the simulation of random heterogeneous medias containing polygonal inclusions and voids); (4) they are more tolerant to large localized deformations. In the past decades, the polygonal
M
elements have been embedded in different numerical methods to solve various problems, e.g., in the FEM for linear elasticity [38] and finite elasticity [35], in the
ED
smoothed FEM for solid mechanics [39], in the extended FEM [40] and NMM [19] for fracture analysis, in the scaled boundary FEM for ealstodynamics [41] and in the
PT
discontinuous Galerkin method involving heat flow [42]. In this paper, to further enlarge the thermal element library and reveal the
CE
properties of the polygonal elements in the NMM background, the Wachspress polygonal elements [43], are further extended to handle 2D transient heat conduction
AC
problems within the platform of the NMM. To this end, the rest of the paper is organized as follows. In Section 2, the governing equations are provided. Section 3 addresses the details about the polygonal NMM for the concerned problems. To verify the proposed method, several numerical examples are tested in Section 4. Finally in Section 5, the concluding remarks are made. 2. Statement of the problem Consider transient heat conduction in a homogeneous isotropic planar domain 4
ACCEPTED MANUSCRIPT
bounded by the contour T
q shown in Fig. 1. The governing partial
differential equation for this kind of problem is [12]
T x, y, t 2T x, y, t 2T x, y, t c k k Q t x 2 y 2
(1)
where ρ, c and k are, respectively, the density, the specific heat at constant pressure
CR IP T
and the thermal conductivity of the material. denotes partial derivative. T, t and Q stand for, respectively, the temperature, the time and the heat source. (x, y) are the Cartesian coordinates plotted in Fig. 1. q
n
AN US
q
y
x
o
M
T
T
ED
Fig. 1. Transient heat conduction in a 2D domain
The associated boundary conditions are
T T nx k ny q on q x y
(3)
CE
k
(2)
PT
T T on T
AC
where T and q are, respectively, the prescribed temperature on Dirichlet boundary
T and the enforced heat flux on Neumann boundary q . (nx , ny ) n is the
outward unit normal to the domain as illustrated in Fig. 1. As a time-dependent problem, the initial condition should also be provided, i.e.,
T x, y,0 T0 ( x, y)
(4)
3. Unsteady heat conduction simulation by the polygonal NMM 3.1 A sketch of the NMM 5
ACCEPTED MANUSCRIPT
The NMM was originally brought out by Shi [14] to study rock engineering problems with satisfactory efficiency and accuracy. The predominance of the NMM, versus other frequently used computational methods (e.g., the FEM and the BEM) results from its dual cover system, i.e., the mathematical cover (MC) and the physical cover (PC). To discretize a physical domain by the NMM, the MC is firstly generated by a series of mathematical patches (MPs) composed of user-defined mathematical
CR IP T
elements. Different from the finite element (FE) mesh, the MC may be inconsistent with the boundary, as long as it is large enough to cover the problem domain. Then, the physical patches (PPs) are obtained by the intersection of MPs and the physical region. All the PPs form the PC. Finally, the manifold elements (MEs) are produced
AN US
through the common zone of PPs.
To further demonstrate the above process, the discretization of the domain
1
2 in Fig. 2a is presented. To cover the whole region, the MC made up of
two hexagonal MPs, that is, M1 and M2 in Fig. 2b, is chosen. The intersection of the
M
physical domain and the MPs (see Fig. 2c) gives 4 PPs in Fig. 2d, that is, P1 and P2 from M1 and , P3 and P4 from M2 and . Above 4 PPs finally generate 6 MEs,
ED
labeled by E1 (P1), E2 (P1, P3), E3 (P3), E4 (P2), E5 (P2, P4) and E6 (P4) in Fig. 2e, where
1
M1
M2
CE
2
PT
the bracketed contents represent the connected PPs.
(b)
AC
(a)
(c)
P1
P3
P2
P4
(d) 6
ACCEPTED MANUSCRIPT
P1
P3
E1
E2
E3
P2
P4
E4
E5
E6
(e) Fig. 2. Illustration of the NMM discretization: (a) the physical domain, (b) the MC, (c) intersection of the physical domain and MPs, (d) generation of the PPs and (e) formations of the MEs
CR IP T
Based on the above concepts, the NMM field interpolation function can be obtained. For the present problem, the temperature in a certain ME E is approximately written as n
T h ( x, y, t ) wi ( x, y )T i ( x, y, t )
(5)
AN US
i 1
where n is the amount of PPs shared by E. wi ( x, y) is the weight function on the ith PP, inherited from that defined on the MP MI containing the very PP, i.e., wI ( x, y) and will be further discussed in the coming subsection. T i ( x, y, t ) is the local
M
function defined on the ith PP and is commonly chosen as
(6)
ED
T i ( x, y, t ) P( x, y)ai ( x, y, t )
where a i is the vector of the generalized degrees of freedom (DOFs) associated with
PT
the ith PP. P( x, y) is the polynomial basis matrix, i.e., P( x, y) 1, x, y, x 2 , xy, y 2 , ...
CE
(7)
3.2 Weight functions for the polygonal elements
AC
According to the NMM fundamentals [14], the weight function wI ( x, y) on the
MP MI should meet the following conditions, that is,
wI ( x, y) C 0 , ( x, y) M I
(8)
wI ( x, y) 0, ( x, y) M I
(9)
where C 0 is the set of all continuous functions.
7
ACCEPTED MANUSCRIPT
Besides, to assure the conforming property, the weight functions on all the MPs containing the point ( x, y) should still satisfy
wJ ( x, y ) 1
(10)
J ( x , y )M J
Although, the shape of mathematical elements may theoretically be user-defined,
CR IP T
to simplify the generation of MC and take advantage of the well-constructed FEs, the noted triangular and quadrilateral elements (or shape functions) are extensively adopted as the NMM mathematical elements (or weight functions) (e.g., see [14, 17, 23]). In the present, the polygonal FEs are further borrowed to solve time-dependent heat transfer problems in the NMM. The birth of the polygonal FEs in the 1970s owes
AN US
to the pioneer, Eugene L. Wachspress and the original intension is to improve the discretization flexibility [43], besides which, more advantages of polygonal elements have been gradually revealed as mentioned in Section 1. In recent years, after the work of Wachspress, several other polygonal FE shape functions have been proposed,
M
including the mean value coordinates by Floater [44], the maximum entropy approximation by Sukumar [45], the natural neighbor shape functions by Sukumar
ED
and Tabarraei [37] and the metric coordinates by Malsch et al. [46]. The comparisons between different polygonal FE shape functions were comprehensively performed by
PT
Sukumar and Malsch [47].
Basically, the above mentioned polygonal shape functions all meet the conditions
CE
listed in Eqs. (8)-(10). In view that it is the earliest and also the simplest form on convex elements [47], the rational Wachspress FE shape function is adopted as the
AC
NMM weight function. For the N-sided convex polygonal element in Fig. 3, the corresponding weight function at point p ( x, y) associated with the vertex vI is
wI ( x, y )
I ( x, y )
(11)
N
J 1
J
( x, y )
where
8
ACCEPTED MANUSCRIPT
I x, y
A vI 1 , vI , vI 1 cot I cot I 2 A vI 1 , vI , p A vI , vI 1 , p pv
(12)
I
cot I cot I
vI p vI vI 1
(13)
vI p vI vI 1 vI vI 1 vI p
(14)
vI vI 1 vI p
CR IP T
where the last expression in Eq. (12) was contributed by Meyer et al. [48].
A p1 , p2 , p3 denotes the area bounded by three vertices p1 , p2 and p3 . I is the included angle between the segments vI p and vI vI 1 , and I is that between
product and module.
AN US
vI vI 1 and vI p . “ ”, “ ” and “| |” are, respectively, the vector dot product, cross
N
1 2
N-1
ED
M
p( x, y)
I I
vI
vI 1
PT
Fig. 3. Illustration of an N-sided polygonal element
What should be further emphasized is that for the case when p lies on certain
CE
element edge, i.e., on side vI vI 1 or vI vI 1 , Eq. (12) is invalid. Under this
AC
circumstances, wI ( x, y) can be calculated through interpolation owing to its linear variation along element edges. 3.3 NMM global equations Presently, the NMM discrete equations for the concerned problems are derived by virtue of the modified variational principle [49]. The equivalent functional in correspondence with Eqs. (1) and (3) is 2 2 T 1 T 1 T (T ) cT k + k QT d qTd q t 2 x 2 y
9
(15)
ACCEPTED MANUSCRIPT
Further, seeing that the MC may not conform to the domain boundaries and the generalized DOFs in Eq. (6) may possess no physical meaning, we utilize the penalty method, adopted widely in the NMM literatures (e.g., see [14, 17, 26, 50]), to approximately impose the Dirichlet boundary condition in Eq. (2). Accordingly, we obtain the modified functional as
1 (T T ) (T T )d 2 T
(16)
CR IP T
* (T ) (T )
where is the penalty.
Based on the variational principle, only the true solution T can satisfy the extreme condition of * , i.e., * (T ) 0 , through which and Eqs. (5), (6) as well, the
AN US
NMM global formulations for transient heat conduction are finally deduced as KT + CT F
(17)
where T and T are, respectively, the vector of DOFs and their time derivatives. K, C and F represent, respectively, the thermal conductivity matrix, the heat capacity
M
matrix and the equivalent thermal load vector. They are firstly computed on individual ME and then summed. For some ME E,
K E E BT kBd
FE
E
(18)
T wi P c wi P d
wi P
T
Qd
PT
E
T wi P wi P d
ED
CE
TE
TE
(19)
T T wi P Td wi P qd E q
(20)
CE
where the superscript T symbolizes the matrix transpose. E is the domain occupied by E. TE and qE denote, respectively, the essential and natural boundary
AC
associated with E. The matrix B can be further expressed as
B B1 B2 ... Bi ... Bn
(21)
( wi P), x Bi ( wi P), y
(22)
with
3.4 Computer implementations
10
ACCEPTED MANUSCRIPT
To solve the linear system of the ordinary differential equations in Eq. (17), the commonly used Euler backward difference scheme [49] is adopted. Accordingly, supposing that the DOFs vector Tt at the instant t is known and the time step size is t , we can obtain Tt t at t t by
C C K t t Tt t Ft t Tt t t
CR IP T
(23)
To acquire K t t , C and Ft t in Eq. (23), the associated spatial integration scheme should also be provided. In the original NMM [14], simplex integration scheme is proposed to analytically compute the domain integration of polynomial
AN US
integrands. In view of that the Wachspress-type weight functions for n-gons (n>4) are rational, this scheme is infeasible herein. In this sense, the numerical integration strategy is used. Recent years, many studies have been carried out with regard to the numerical integration on polygonal elements, e.g., see [51-54]. In the present paper, a
M
conventional technique adopted in [40, 55] is taken. In the computation, every ME is firstly triangulated, and then regular Gaussian quadrature rules are applied on each
ED
triangle. Finally, the results on all the sub-elements add up to those on the associated ME.
PT
4. Numerical examples
In this section, to verify the feasibility and accuracy of the proposed method, four
CE
transient heat conduction cases, under Dirichlet boundary condition or mixed Dirichlet and Neumann boundary conditions, are conducted. As for the shape of the
AC
MP, theoretically, any arbitrarily shaped mathematical element is available to form the MPs, however, in light of that the MC can be inconsistent with the physical boundary and the accuracy of regular elements [56], the equilateral triangles and squares are mostly used in the NMM, e.g., as in [14, 17, 18, 23, 28]. In this work, on the ground that among all the regular n-gons (n>4), only hexagons can individually cover the plane, MPs composed of regular hexagonal mathematical elements are adopted throughout the analysis. Furthermore, to clearly quantify the discretization, besides the amount of PPs and MEs, the element size (the side length of a mathematical element) is also used. In addition, to avoid the rank deficiency trouble in the use of 11
ACCEPTED MANUSCRIPT higher-order local function [57, 58], the polynomial basis in Eq. (7) is chosen to be constant. 4.1 A square plate under Dirichlet boundary condition In the first example, heat conduction in a square plate without heat source is considered. As shown in Fig. 4, the width of the plate is L. The temperature T is enforced along the entire boundary and the initial condition is set such that T ( x, y,0) 10sin( x)sin( y) .
CR IP T
In the analysis, the following parameters are taken: L , 1.0, c 1.0 k 1.0 and T 0 . Accordingly, the theoretical temperature field is [12]
T ( x, y, t ) 10sin( x)sin( y)exp(2t )
(24)
the Fourier’s law as qx 10cos( x)sin( y ) exp(2t ) q y 10sin( x) cos( y ) exp(2t ) y
ρ, c, k
T
T
L
x
L
(a)
(b)
(c)
(d)
AC
CE
o
PT
A
ED
B
T
(25)
M
T
AN US
and the components of the heat flux vector, that is, qx and qy, can be derived through
Fig. 4. Transient heat conduction in a square plate under Dirichlet boundary condition: (a) physical model, (b) discretization when h=0.1, (c) discretization when h=0.2 and (d) discretization when h=0.4 12
ACCEPTED MANUSCRIPT
Firstly, the sensitivity study on the effect of the penalty number is investigated. For this purpose, the variations of temperature along the line y =π/2 are concerned with ∆t = 0.02, under the discretization in Fig. 4b (element size h = 0.1, 959 PPs and 429 MEs) at different λ, that is, λ =10k, 102k, 103k，104k，106k，108k and 1010k. The simulated temperatures when x =iπ/12 (i= 1, 2…6) at two instants, i.e., t =0.5 and 1.0, are, respectively, shown in Table 1 and 2, together with the exact solutions from Eq.
CR IP T
(24). Accordingly, we can find that the error of the present results with relatively small λ (e.g., 10k and 102k) is greater than that of larger λ (e.g., 106k and more) at a given point; besides, for the cases when λ>104k, the difference among NMM solutions is ignorable. On this ground, λ will be fixed at 1010k in the following computations.
λ/k =10
1.017469
1.909451
λ/k =102
0.970114
1.866591
3
0.965342
1.862248
λ/k =104
0.964864
1.861813
λ/k =106
0.964812
1.861766
8
0.964811
λ/k =1010
0.964811
Exact
0.952144
x=5π/12
x=π/2
2.676400
3.270243
3.634246
3.757009
2.638406
3.236455
3.603165
3.726853
2.634538
3.233001
3.599979
3.723760
2.634150
3.232655
3.599660
3.723450
2.634108
3.232617
3.599625
3.723416
1.861765
2.634107
3.232616
3.599624
3.723415
1.861765
2.634107
3.232616
3.599624
3.723415
1.839401
2.601305
3.185934
3.553445
3.678794
PT
λ/k =10
x=π/3
M
x=π/6
λ/k =10
x=π/4
ED
x=π/12
AN US
Table 1 Calculated temperature at selected points under different penalty values when t=0.5
Table 2 Calculated temperature at selected points under different penalty values when t=1.0 x=π/6
x=π/4
x=π/3
x=5π/12
x=π/2
λ/k =10
0.388076
0.728180
1.020464
1.246647
1.385214
1.431934
λ/k =102
0.363929
0.700224
0.989741
1.214063
1.351606
1.397997
3
0.361531
0.697431
0.986660
1.210788
1.348223
1.394579
λ/k =104
0.361291
0.697152
0.986352
1.210460
1.347884
1.394237
λ/k =106
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
8
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
10
0.361265
0.697121
0.986318
1.210424
1.347847
1.394199
0.350274
0.676678
0.956967
1.172040
1.307240
1.353353
CE
x=π/12
AC
λ/k =10
λ/k =10
λ/k =10 Exact
Then, the convergence study is performed concerning the effect of element size. Another two MCs with h=0.2 and 0.4 (see the discretization in Fig. 4c with 262 PPs 13
ACCEPTED MANUSCRIPT
and 110 MEs when h=0.2, and Fig. 4d with 89 PPs and 33 MEs when h=0.4) are further used under ∆t =0.02. The calculated temperatures at two sample points in Fig. 4a, that is, point A: (π/4, π/4) and point B: (π/2, π/2), are respectively, summarized in Fig. 5a and 5b, together with the exact ones. From these figures, we can find that with the decrease of h, the NMM results gradually converge to the analytical solutions. 10.0
Exact NMM-h=0.4 NMM-h=0.2 NMM-h=0.1
8.0 7.0
Temperature
Temperature
4.0
Exact NMM-h=0.4 NMM-h=0.2 NMM-h=0.1
9.0
CR IP T
5.0
3.0
2.0
6.0 5.0 4.0 3.0
1.0
2.0 1.0
0.0 1.0
2.0
3.0
4.0
AN US
0.0
0.0
5.0
Time
0.0
1.0
2.0
3.0
4.0
5.0
Time
(a) (b) Fig. 5. Computed temperatures of sample points at different instants and element size: (a) point A and (b) point B
Next, the effect of time step size on the solution accuracy is investigated. To this
M
end, three time steps, i.e., ∆t =0.1, 0.05 and 0.02, are sequentially examined. The
ED
discretization in Fig. 4b is used. The computed temperatures of point A and B at different instants are plotted in Fig. 6, from which the decay of temperature with time
PT
can be detected. For comparison, the exact results from Eq. (24) are also illustrated therein. Obviously, the temperatures at all time steps match well with the analytical
CE
ones; what’s more, with the decrease of time step, our results are getting closer to the exact values, which also conforms to the convergence rule of the backward difference
AC
scheme.
10.0
5.0
Exact NMM-time step= 0.10 NMM-time step= 0.05 NMM-time step= 0.02
8.0 7.0
Temperature
Temperature
4.0
3.0
2.0
Exact NMM-time step=0.10 NMM-time step=0.05 NMM-time step=0.02
9.0
6.0 5.0 4.0 3.0
1.0
2.0 1.0
0.0
0.0 0.0
1.0
2.0
3.0
4.0
0.0
5.0
1.0
2.0
3.0
Time
Time
14
4.0
5.0
ACCEPTED MANUSCRIPT (a) (b) Fig. 6. Computed temperatures of sample points at different instants and time steps: (a) point A and (b) point B
Lastly, the accuracy of heat flux is further tested under the discretization in Fig. 4b and ∆t = 0.02. The flux components at point A are plotted in Fig. 7, together with the exact ones by Eq. (25). Upon these figures, we can conclude that the present flux
CR IP T
solutions also conform well to the analytical values. -5.0
-5.0
Exact NMM
Exact NMM
-3.0
-3.0
qx
qy
-4.0
-4.0
-2.0
-1.0
-1.0
AN US
-2.0
0.0
0.0 0.0
1.0
2.0
3.0
4.0
0.0
5.0
Time
1.0
2.0
3.0
4.0
5.0
Time
Fig. 7. Computed heat flux of point A at different instants: (a) qx and (b) qy
4.2 A square plate under mixed boundary conditions
M
In this example, we consider a square plate of width L with zero initial
ED
temperature in Fig. 8a. Just as in Example 1, the heat source is also ignored. The temperature T is imposed along the left boundary and the heat flux q is prescribed
PT
on the bottom boundary with the rest boundaries insulated. Note that this problem was previously solved by Li et al. with the meshless local Petrov-Galerkin method (MLPG)
CE
[12].
AC
y
T
o
ρ, c, k
B
q
A
L
x
L
15
ACCEPTED MANUSCRIPT (a)
(b)
Fig. 8. Transient heat conduction in a square plate under mixed boundary conditions: (a) physical model and (b) the discretization when h=2.8
In the modeling, the associated parameters are chosen as: L 100, 1.0,
c 1.0, k 1000, T 0 and q 1000 . Mathematical element of size h=2.8 is applied to form the MPs and the discretization finally generates 1158 PPs and 525 MEs as displayed in Fig. 8b. For comparison, the same time step as in [12], i.e.,
CR IP T
t 0.05 is used. Two sample points with their coordinates: (100, 0) and (100, 50)
(denoted, respectively, as point A and B in Fig. 8a), are concerned. The computed temperatures when t [0, 2.0] for both points are, respectively, depicted in Fig. 9a and 9b, together with those under regular nodes in [12] (note that the MLPG results
AN US
were approximately extracted from Figs. 9 and 10 of reference [12] ). It is easy to achieve that our solutions agree well with the existing ones. 16.0
50.0
14.0
45.0 40.0
12.0
MLPG Present
10.0
25.0 20.0
10.0 5.0 0.0 -5.0 -0.2
0.0
0.2
0.4
0.6
0.8
ED
15.0
1.0
1.2
1.4
1.6
1.8
Temperature
30.0
M
Temperature
35.0
MLPG Present
8.0 6.0 4.0 2.0 0.0
2.0
-2.0 -0.2
2.2
0.0
0.2
0.4
0.6
0.8
Time
1.0
1.2
1.4
1.6
1.8
2.0
2.2
PT
Time
(a)
(b)
CE
Fig. 9. Computed temperatures of sample points at different instants of: (a) point A and (b) point B
4.3 A concave domain under time-dependent boundary condition
AC
Herein, a concave domain with heat source subjected to time-dependent essential
boundary condition is considered. The dimensions and material properties are denoted in Fig. 10a. The temperature T T (t ) is applied on the left boundary and other boundaries are insulated. As for the initial condition, T ( x, y,0) 0 .
16
ACCEPTED MANUSCRIPT y
5 4
ρ, c, k
3
7
Q
12 11
6
2
10
1
9
H
CR IP T
T (t )
13
8
x
W2
W1
AN US
(a)
W1
M
(b)
ED
Fig. 10. Transient heat conduction in a concave domain under time-dependent boundary condition: (a) physical model and (b) the discretization when h=0.02
In the simulation, W1 0.3 , W2 0.4 , H 0.6 , 1.0 , c 1.0, k 1.0 , Q 10 ,
PT
T (t ) 60t , t 0.02 and h =0.02. The discretized domain is given in Fig. 10b, which
involves 1280 PPs and 570 MEs. Temperatures at 13 sample points (see Table 3 for
CE
their coordinates) in Fig. 10a are tested and the computed results at t =1.0 and 10.0 are summarized in Table 3. Concerning that the same problems were tackled by Yao et al.
AC
with the precise integration boundary element method (PIBEM) [59], we also report their solutions in this table. As we can see, the present results are in nice consistence with the reference values and the maximum relative error (RE) is less than 0.53% for t =1.0 and 0.07% for t =10.0.
17
ACCEPTED MANUSCRIPT
Table 3 Calculated temperature at sample points by NMM and PIBEM temperature
sample point t =1.0
t =10.0
NMM
PIBEM
RE (%)
NMM
PIBEM
RE (%)
1: (0.15, 0.1)
55.7389
55.6649
0.133
595.3105
595.2820
0.005
2: (0.15, 0.2)
55.0830
55.0101
0.132
594.5542
594.5438
0.002
3: (0.15, 0.3)
54.4254
54.3731
0.096
593.7956
593.8255
0.005
4: (0.15, 0.4)
54.0200
53.9746
0.084
593.3275
593.3760
0.008
5: (0.15, 0.5)
53.8201
53.7779
0.078
593.0964
593.1536
0.010
6: (0.50, 0.3)
42.4838
42.3573
0.299
580.0919
580.2965
0.035
7: (0.50, 0.4)
42.4340
42.2967
0.325
580.0375
580.2317
0.033
8: (0.50, 0.5)
42.3966
42.2550
0.335
579.9968
580.1868
0.033
9: (0.85, 0.1)
34.6041
34.4424
0.469
570.8062
571.1672
0.063
10: (0.85, 0.2)
35.0822
34.9159
0.476
571.3773
571.7179
0.060
11: (0.85, 0.3)
35.5935
35.4115
0.514
571.9888
572.2947
0.053
12: (0.85, 0.4)
35.9529
35.7641
0.528
572.4200
572.7055
0.050
13: (0.85, 0.5)
36.1552
35.9660
0.526
572.6634
572.9410
0.048
ED
AN US
CR IP T
(x, y)
M
number:
4.4 A square domain with multiple holes
PT
In the last example, to further demonstrate the advantages of the proposed method in discretization and accuracy as well, a square domain of width W with five holes in
CE
Fig. 11a is studied. The side length of the diamond (square) hole is a while the circular holes are of the same radius b. Material parameters and other dimensions are
AC
all presented therein. The boundary conditions are prescribed such that the temperature along the left boundary is T , the heat flux q1 (t ) and q2 (t ) continuously enter the plate from the top side and the central hole boundaries, respectively, and the remaining edges are all insulated.
18
ACCEPTED MANUSCRIPT q1 (t )
ρ, c, k
W1
B
b
b y
o
C
T
q2
q2
a
b
W2
x
A
q2
W2
b
D
W1
Q(t) W1
W2
W2
W1
W
(a)
CR IP T
q2
(b)
AN US
Fig. 11. Transient heat conduction in a square domain with five holes: (a) physical model and (b) the MC when h=0.15
In the NMM computations, the corresponding parameters are chosen as:
W 10.0, W1 W2 2.5 , a 1.5, b 1.0 , 1.0 , c 1.0,
k 1.0 , T 0 ,
q1 10(1 t ) and q2 5 . At the same time, the initial temperatures are set to be zero
M
and the heat source is Q(t ) 10sin(2t ) . Regular MC with h =0.15 as in Fig. 11b is
ED
used, regardless the size, location and boundaries of all interior holes and the exterior domain boundary as well. After some geometrical operations, the final NMM
AC
CE
PT
discretization is obtained in Fig. 12a, which includes 3457 PPs and 1567 MEs.
(a)
(b)
Fig. 12. Discretization by NMM and FEM: (a) NMM discretization and (b) FEM mesh
19
ACCEPTED MANUSCRIPT
To the authors’ knowledge, there are no reference solutions to this problem. Therefore, our results are compared with those from the widely used FE software ANSYS (version R14.5). The corresponding FE mesh is shown in Fig. 11b, with 2744 nodes and 2476 4-noded quadrilateral elements. From Figs. 11b, 12a and 12b, the unique feature of the NMM in meshing, as addressed in Section 1 and 3, can be easily inferred. In the comparative study, the time step in both methods is fixed at t 0.02 .
CR IP T
In addition, the Newmark scheme with ANSYS default parameters is used in the FEM analysis. The temperatures at 4 sample points in Fig. 11a, i.e., point A: (2.5, 0.0), point B: (0.0, 2.5), point C: (-2.5, 0.0) and point D: (0.0, -2.5), are, respectively, examined. The calculated results are plotted in Fig. 13a~d. As we can see, the
AN US
difference of temperatures by the NMM and the FEM is almost indistinguishable at all points.
45
18
FEM NMM
16
35
14 12 10 8 6
2 0
0.0
PT
-2 1.0
2.0
3.0
4.0
25 20 15 10
ED
4
Temperature
30
M
Temperature
FEM NMM
40
5 0 -5 5.0
0.0
1.0
2.0
Time
(a)
CE
16
4.0
5.0
(b) 16
FEM NMM
14
FEM NMM
14
12
12
10
10
Temperature
AC Temperature
3.0
Time
8 6 4
8 6 4
2
2
0
0
-2 0.0
1.0
2.0
3.0
4.0
-2
5.0
0.0
Time
1.0
2.0
3.0
Time
(c)
(d)
20
4.0
5.0
ACCEPTED MANUSCRIPT Fig. 13. The calculated temperatures of sample points at different instants by NMM and FEM: (a) point A, (b) point B, (c) point C and (d) point D
5. Concluding remarks In this work, the numerical manifold method was developed to tackle 2D unsteady heat transfer problems on hexagonal elements. Besides the description of governing equations, the introduction of the NMM, and the derivation of discrete formulations,
CR IP T
details about the weight functions of the Wachspress polygonal elements, the solution strategy to the time-dependent system of equations and the numerical integration schemes were also provided. Four numerical examples with increasing complexity are tested, under Dirichlet boundary condition or mixed Dirichlet and Neumann boundary conditions. Mathematical covers composed of regular hexagonal mathematical
AN US
elements were adopted throughout the simulation. The effect of some key parameters, i.e., the penalty value, the element size and the time step, on the results was investigated. The computed thermal fields agree well with the available reference solutions.
M
For the considered continuous heat transfer problems, compared with some other represented numerical tools, e.g., the FEM, the BEM and the meshless methods, the
ED
dominant advantage of the proposed approach lies in discretization, that is, the mathematical cover can be inconsistent with both external and internal boundaries
PT
(e.g., see example 4), which may reduce the meshing cost in part and exploit the accuracy of regular elements as well. Although only regular hexagons were used in
CE
this study, all convex n-gons (n>4) and their combinations are welcome from the essence of the Wachspress type polygonal NMM. The work on arbitrary polygonal
AC
elements is proceeding and the results are awaited. Acknowledgements The present work was supported by the National Natural Science Foundation of China (Grant Nos. 11462014 and 11572282), the Provincial Natural Science Foundation of Jiangxi, China (Grant No. 20151BAB202003), the Science and Technology Program of Educational Committee of Jiangxi Province of China (Grant Nos. GJJ14526 and GJJ150752). 21
ACCEPTED MANUSCRIPT Reference [1] M.N. Ozisik, Boundary values problems of heat conduction, Internatinal Textbook Company, Scranton, 1968. [2] S. Singh, P.K. Jain, Rizwan-uddin, Analytical solution to transient heat conduction in polar coordinates with multiple layers in radial direction, Int J Therm Sci, 47 (2008) 261-273. [3] H.J. Jiang, H.L. Dai, Analytical solutions for three-dimensional steady and transient heat conduction problems of a double-layer plate with a local heat source, Int J Heat Mass Tran, 89 (2015) 652-666. [4] W.Z. Dai, R. Nassar, A finite difference scheme for solving the heat transport equation at the
CR IP T
microscale, Numer Meth Part D E, 15 (1999) 697-708.
[5] C.C. Wang, Application of the maximum principle for differential equations in combination with the finite difference method to find transient approximate solutions of heat equations and error analysis, Numer Heat Tr B-Fund, 55 (2009) 56-72.
[6] T.D. Hien, M. Kleiber, Stochastic finite element modelling in linear transient heat transfer, Comput Method Appl M, 144 (1997) 111-124.
AN US
[7] Y.X. Chen, X.B. Liang, Y. Liu, J.Y. Bai, B.S. Xu, Finite element modeling of coating formation and transient heat transfer in the electric arc spray process, Int J Heat Mass Tran, 53 (2010) 2012-2021. [8] S.C. Mishra, H.K. Roy, Solving transient conduction and radiation heat transfer problems using the lattice Boltzmann method and the finite volume method, J Comput Phys, 223 (2007) 89-107. [9] M.J. Guillot, S.C. McCool, Effect of boundary condition approximation on convergence and accuracy of a finite volume discretization of the transient heat conduction equation, Int J Numer Method H, 25 (2015) 950-972.
M
[10] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method, Eng Anal
ED
Bound Elem, 26 (2002) 119-132.
[11] A.I. Abreu, A. Canelas, W.J. Mansur, A CQM-based BEM for transient heat conduction problems in homogeneous materials and FGMs, Appl. Math. Model., 37 (2013) 776-792.
PT
[12] Q.H. Li, S.S. Chen, G.X. Kou, Transient heat conduction analysis using the MLPG method and modified precise time step integration method, J Comput Phys, 230 (2011) 2736-2750. [13] T. Sophy, A. Da Silva, A. Kribeche, H. Sadat, An alternative space-time meshless method for
CE
solving transient heat transfer problems with high discontinuous moving sources, Numer Heat Tr B-Fund, 69 (2016) 377-388. [14] G.H. Shi, Manifold method of material analysis,
Transcations of the 9th Army Confernece on
AC
Applied Mathematics and ComputingMinneapolis, Minnesota, 1991, pp. 57-76. [15] G.W. Ma, X.M. An, L. He, The numerical manifold method: a review, Int J Comp Meth-Sing, 7 (2010) 1-32. [16] R.J. Tsay, Y.J. Chiou, W.L. Chuang, Crack growth prediction by manifold method, J Eng Mech-Asce, 125 (1999) 884-890. [17] G.W. Ma, X.M. An, H.H. Zhang, L.X. Li, Modeling complex crack problems using the numerical manifold method, International Journal of Fracture, 156 (2009) 21-35. [18] H.H. Zhang, L.X. Li, X.M. An, G.W. Ma, Numerical analysis of 2-D crack propagation problems using the numerical manifold method, Eng Anal Bound Elem, 34 (2010) 41-50. [19] H.H. Zhang, S.Q. Zhang, Extract of stress intensity factors on honeycomb elements by the numerical manifold method, Finite Elem Anal Des, 59 (2012) 55-65. 22
ACCEPTED MANUSCRIPT
[20] X.M. An, Z.Y. Zhao, H.H. Zhang, L. He, Modeling bimaterial interface cracks using the numerical manifold method, Eng Anal Bound Elem, 37 (2013) 464-474. [21] H. Zheng, D.D. Xu, New strategies for some issues of numerical manifold method in simulation of crack propagation, Int J Numer Meth Eng, 97 (2014) 986-1010. [22] H.H. Zhang, G.W. Ma, F. Ren, Implementation of the numerical manifold method for thermo-mechanical fracture of planar solids, Eng Anal Bound Elem, 44 (2014) 45-54. [23] H.H. Zhang, G.W. Ma, Fracture modeling of isotropic functionally graded materials by the numerical manifold method, Eng Anal Bound Elem, 38 (2014) 61-71. [24] J. He, Q.S. Liu, G.W. Ma, B. Zeng, An improved numerical manifold method incorporating hybrid
CR IP T
crack element for crack propagation simulation, International Journal of Fracture, 199 (2016) 21-38.
[25] H. Zheng, F. Liu, X.L. Du, Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method, Comput Method Appl M, 295 (2015) 150-171.
[26] G.X. Zhang, Y. Zhao, X.C. Peng, Simulation of toppling failure of rock slope by numerical manifold method, Int J Comp Meth-Sing, 7 (2010) 167-189.
[27] Y.J. Ning, X.M. An, G.W. Ma, Footwall slope stability analysis with the numerical manifold
AN US
method, Int J Rock Mech Min, 48 (2011) 964-975.
[28] L.N.Y. Wong, Z.J. Wu, Application of the numerical manifold method to model progressive failure in rock slopes, Eng Fract Mech, 119 (2014) 1-20.
[29] L.F. Fan, X.W. Yi, G.W. Ma, Numerical manifold method (Nmm) simulation of stress wave propagation through fractured rock mass, Int J Appl Mech, 5 (2013) 238-249.
[30] G.F. Zhao, X.B. Zhao, J.B. Zhu, Application of the numerical manifold method for stress wave propagation across rock masses, Int J Numer Anal Met, 38 (2014) 92-110.
M
[31] X.L. Qu, Y. Wang, G.Y. Fu, G.W. Ma, Efficiency and accuracy verification of the explicit numerical manifold method for dynamic problems, Rock Mech Rock Eng, 48 (2015) 1131-1142. [32] Q.H. Jiang, S.S. Deng, C.B. Zhou, W.B. Lu, Modeling unconfined seepage flow using
ED
three-dimensional numerical manifold method, J Hydrodyn, 22 (2010) 554-561. [33] H. Zheng, F. Liu, C.G. Li, Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method, Appl. Math. Model., 39 (2015) 794-808.
PT
[34] Y. Wang, M.S. Hu, Q.L. Zhou, J. Rutqvist, A new second-order numerical manifold method model with an efficient scheme for analyzing free surface flow with inner drains, Appl. Math. Model., 40
CE
(2016) 1427-1445.
[35] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, Polygonal finite elements for finite elasticity, Int J Numer Meth Eng, 101 (2015) 305-328.
AC
[36] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, Comput Method Appl M, 306 (2016) 216-251.
[37] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Int J Numer Meth Eng, 61 (2004) 2045-2066.
[38] A. Tabarraei, N. Sukumar, Application of polygonal finite elements in linear elasticity, Int J Comp Meth-Sing, 3 (2006) 503-520. [39] K.Y. Dai, G.R. Liu, T.T. Nguyen, An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics, Finite Elem Anal Des, 43 (2007) 847-860. [40] A. Tabarraei, N. Sukumar, Extended finite element method on polygonal and quadtree meshes, Comput Method Appl M, 197 (2008) 425-438. 23
ACCEPTED MANUSCRIPT
[41] Z.H. Zhang, Z.J. Yang, J.H. Li, An adaptive polygonal scaled boundary finite element method for elastodynamics, Int J Comp Meth-Sing, 13 (2016). [42] J. Jaskowiec, P. Plucinski, A. Stankiewicz, Discontinuous Galerkin method with arbitrary polygonal finite elements, Finite Elem Anal Des, 120 (2016) 1-17. [43] E.L. Wachspress, A rational finite element basis, Academic Press, New York, 1975. [44] M.S. Floater, Mean value coordinates, Comput Aided Geom D, 20 (2003) 19-27. [45] N. Sukumar, Construction of polygonal interpolants: a maximum entropy approach, Int J Numer Meth Eng, 61 (2004) 2159-2181. polygons, Journal of Graphics Tools, 10 (2005) 27-39.
CR IP T
[46] E.A. Malsch, J.J. Lin, G. Dasgupta, Smooth two-dimensional interpolations: a recipe for all [47] N. Sukumar, E.A. Malsch, Recent advances in the construction of polygonal finite element interpolants, Arch Comput Method E, 13 (2006) 129-163.
[48] M. Meyer, A. Barr, H. Lee, M. Desbrun, Generalized barycentric coordinates on irregular polygon, Journal of Graphics Tools, 7 (2002) 13-22.
[49] J.F. Wang, Y.M. Cheng, A new complex variable meshless method for transient heat conduction
AN US
problems, Chinese Phys B, 21 (2012).
[50] H.H. Zhang, G.W. Ma, L.F. Fan, Thermal shock analysis of 2D cracked solids using the numerical manifold method and precise time integration, Eng Anal Bound Elem, 75 (2017) 46-56. [51] G. Dasgupta, Integration within polygonal finite elements, J Aerospace Eng, 16 (2003) 9-18. [52] S. Natarajan, S. Bordas, D.R. Mahapatra, Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping, Int J Numer Meth Eng, 80 (2009) 103-134. [53] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on
M
irregular convex polygons and polyhedrons, Comput Mech, 47 (2011) 535-554. [54] E.B. Chin, J.B. Lasserre, N. Sukumar, Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra, Comput Mech, 56 (2015) 967-981.
ED
[55] A.R. Khoei, R. Yasbolaghi, S.O.R. Biabanaki, A polygonal finite element method for modeling crack propagation with minimum remeshing, International Journal of Fracture, 194 (2015) 123-148. [56] H.H. Zhang, Y.L. Chen, L.X. Li, X.M. An, G.W. Ma, Accuracy comparison of rectangular and
PT
triangular mathematical elements in the numerical manifold method, Analysis of Discontinuous Deformation: New Developments and Applications, (2010) 297-303.
CE
[57] R. Tian, G. Yagawa, H. Terasaka, Linear dependence problems of partition of unity-based generalized FEMs, Comput Method Appl M, 195 (2006) 4768-4782. [58] X.M. An, L.X. Li, G.W. Ma, H.H. Zhang, Prediction of rank deficiency in partition of unity-based
AC
methods with plane triangular or quadrilateral meshes, Comput Method Appl M, 200 (2011) 665-674. [59] W.A. Yao, B. Yu, X.W. Gao, Q. Gao, A precise integration boundary element method for solving transient heat conduction problems, Int J Heat Mass Tran, 78 (2014) 883-891.
24