- Email: [email protected]
hyperbolic heat conduction type problems
International Journal of Heat and Mass Transfer 101 (2016) 365–372
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A consistent Moving Particle System Simulation method: Applications to parabolic/hyperbolic heat conduction type problems Tao Xue a,b, Kumar K. Tamma b,⇑, Xiaobing Zhang a a b
School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing, Jiangsu 210094, PR China Department of Mechanical Engineering, University of Minnesota-Twin Cities, 111 Church St. SE, Minneapolis, MN 55455, United States
a r t i c l e
i n f o
Article history: Received 8 March 2016 Received in revised form 5 May 2016 Accepted 6 May 2016 Available online 1 June 2016 Keywords: MPS Particle based method Parabolic/hyperbolic heat conduction GS4 i-Integration
a b s t r a c t In this paper, a consistent MPS (Moving Particle System Semi-implicit/Moving Particle System Simulation) method is established emanating from Taylor Expansion and Gauss’s Theorem. The proposed MPS method preserves the consistency between the Gradient and Laplacian approximations, and also has the ability to recover the conventional MPS method which lacks consistency. In addition, the proposed method preserves the first order accuracy in space whereas the traditional MPS method loses the accuracy in the case with arbitrary and non-uniform particle discretization. The applicability of the proposed method for the analysis of first- and second-order transient hyperbolic and parabolic systems is illustrated via a recently developed C–F model (Cattaneo–Fourier model) of heat conduction that has been proven to govern thermal transport in different scales, and also capture the various heat conduction physical phenomena. The approaches of imposing different boundary conditions in the heat conduction problems are investigated comprehensively as well. In addition, a unified GS4 i-Integration framework is exploited in the numerical simulations involving first- and second-order systems in time. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The MPS (Moving Particle System Semi-implicit/Moving Particle System Simulation) method is a particle based method under the description of Lagrangian meshless method originally proposed for simulations of incompressible fluid dynamics with the objective to capture the free surface [1]. The convenience of MPS method as well as the other particle based meshless methods such as Smooth Particle Hydrodynamics (SPH) method [2,3], is to provide approximations to the strong form of the partial differential equations (PDEs) via using Gradient operator and Laplacian operator. Nevertheless, the (original) MPS method offers the gradient of a scalar value which is mainly used for calculation of the pressure gradient in the Navier–Stokes equation. On the other hand, the Laplacian operator is designed with regards to a vector value which is fluid velocity; consequently this is not consistent with the gradient operator from the perspective of Gauss’s Theorem and basic differential calculus principles. During the past several decades, the MPS method has been widely implemented in various engineering fields. Koshizuka [1,4] established the basic formulations of the MPS method and applied ⇑ Corresponding author. E-mail addresses: [email protected] (T. (K.K. Tamma), [email protected] (X. Zhang).
Xue),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.05.020 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
[email protected]
it to dam break problems for better capturing the free surface. The enhanced MPS method has been developed to improve the stability and accuracy of the MPS methods [5]. One major challenge of the MPS-type simulation in solving PDEs is the boundary condition implementations, which also exists in the other particle based methods. The issue of boundaries in the MPS applications of fluid dynamics is resolved by adding several layers of wall particles, and these wall particles are also used to prevent the fluid particles from penetrating into the wall. However, the approach of imposing boundary conditions, including Dirichlet boundary conditions and Neumann boundary conditions, has not been fully investigated to-date in the MPS method from the perspective of solving a general PDE. Heat Conduction The theoretical basis of heat conduction and transport relies heavily on the heat flux constitutive models and the fusion of the mechanisms of heat transport with regard to heat conduction. Various heat flux models have been proposed for decades and different heat flux models lead to different types of differential equation governing the heat conduction phenomenon in the solid state [6–8]. To better govern the physics of the heat conduction processes and predict the material properties, a novel C-(Cattaneo-type) and F-(Fourier-type) processes conduction constitutive [9,10] has been developed from the physics of the Boltzmann Transport Equation; and the model acknowledges the notion of the simultaneous coexistence of both slow Cattaneo-type C-processes and fast Fourier-type F-processes in the process of
366
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
heat conduction. Different types of heat flux models, lead to the different heat conduction governing equations and represent different mathematical and physical characteristics. In the C–F model, different heat conduction model numbers F T yield different types of PDEs, including the hyperbolic (F T ¼ 0) and parabolic (F T ¼ 1) heat conduction type problems. Hence, the heat conduction C–F model provides us with a comprehensive mathematical example to test and verify the MPS formulation. The Gradient and Laplacian operators can be directly approximated by the MPS method with particle discretization, and the particles carry the physical entity of temperature and transfer it to the neighboring particles. There have been only a limited or scant literature for solving the heat conduction and transport problems via using the MPS method, and the problems involving different types of boundary effects have not been thoroughly investigated. Therefore, the C–F heat conduction model with different type of boundary conditions offers a numerical environment for improving and testing the MPS method, especially from the aspect of the implementation of the boundaries in the PDEs approximated by the MPS method. Time Integration The particle based methods require neighboring particle detecting process for each particle and it turns out to be the very expensive in practical simulations, especially for large scale systems involving large number of particles. Under this general scenario of computational efficiency, the use of explicit time integration schemes with a minimum of second order accuracy in time are recommended in contrast to the first order explicit Euler-based methods. A unified time integration framework, under the umbrella of the so-called generalized single-step single-solve (GSSSS) family of algorithms, including not only implicit, but also explicit algorithms (the representations can be also written in a linear three-step form), are proposed in [11,12]. This simulation framework has been originally developed in the field of computational structural dynamics. It encompasses most existing second order time accurate algorithms, such as Newmark, Central Difference method and HH-a, and more novel and optimal algorithms and advances. Additionally the implementation to first order systems such as heat conduction problem has been developed and proved to retain the same advantages via an isochronous integration [i-Integration] [13] process within the GSSSS framework. In order to improve and extend the traditional MPS method, the main objectives of this paper focus upon the following: (1) To enhance the MPS method to heat conduction problems, an improved consistent MPS-based method is presented depending upon the established MPS Gradient and Laplacian operators. Via the C–F heat flux constitutive model, the consistent MPS-based method is tested and implemented to illustrate the various heat conduction mechanisms in both the parabolic and hyperbolic differential equations. (2) To investigate the implementation of boundaries in the consistent MPS method, as well as other particle-based method, the approach of imposing different types of boundary effects is presented in detail, and these include imposed temperature (Type 1), heat flux (Type 2), convection (Type 3) and the like. (3) To improve the computational efficiency and the numerical compatibility in terms of the developed MPS based method, the second order time accurate unified GS-4 i-Integration framework is implemented within the numerical simulations with regard to both the first- and second- order transient systems.
of the traditional MPS method. Specifically, the terms related to the Gradient and Laplacian of temperature in heat conduction will be developed from Taylor Series, and the consistency between these two operators is strictly preserved via using Gauss’s Theorem; this is in contrast to past practices. 2.1. Standard MPS method In the MPS method the gradient and the Laplacian operators are formulated as:
hr/ii ¼ D
E
r2 /
i
Ndim X /j /i rj ri W ij n0 j–i jjrj r i jj2 ¼
2Ndim X ð/ /i ÞW ij n0 k0 j–i j
ð1Þ
where / is a physical quantity, Ndim is dimensional number of the simulating space, r is coordinate vector, the subscripts i and j are the targeted particle and its neighboring particles, W ij is the weighting function,1 and n0 denotes the particle number density. The formulations of W ij ; n0 and k0 are listed as follows:
W ij ¼
R=rij 1; 0 6 r ij < R
ð2Þ
R6r
0;
where r ij ¼ r j r i ; R is the radius of the influence domain, and
n0 ¼
X j–i
P W 0ij ;
k0 ¼
0 0 j–i W ij r j
P
2 r0i
0 j–i W ij
ð3Þ
It is worth noting that the physical quantity / is a scalar value in the Gradient operator and a vector value in the Laplacian operator; specifically, the objectives of these two operators are to formulate pressure and velocity in the Navier–Stokes equation which the original MPS was developed for. Therefore, these operators are not consistent with each other and the relationship does not satisfy the divergence theorem. In addition, n0 and k0 are obtained by the initial weighting function, W 0ij . 2.2. Consistent MPS-based method We herein present a consistent formulation. The MPS method discretizes the domain as a collection of particles in such a way the physical quantities at each of the particles are affected by its neighboring particles. Consider a domain in Euclidean space, E. Let r i 2 X0 E denotes the geometric position vector of a certain position in the domain and r j 2 X0 E denotes another point which is close to position r i . Therefore, any physical quantity at each point of space and time can be defined as uðri ; tÞ : X0 I ! E. To simplify the formulations, the physical quantity uðr i ; tÞ is replaced by ui in the rest of this paper, the subscript i represents particle label and the subscript, ij, represents the relationship between particle i and j. Gradient Operator. In mathematics, a Taylor series is a representation of a function as an infinite sum of terms that are calculated from the values of the function’s derivatives at a single point. As mentioned above, the first order Taylor Series of a general physical quantity uj about r i at a fixed time t yields
uj ¼ ui þ ruij ðrj ri Þ þ Oðkrj ri k2 Þ
ð4Þ
In the view of Eq. (4), omitting the higher order terms, we have
2. Theory This section recapitulates the main idea used for the derivations of the Gradient and Laplacian operators of a field property in terms
1 The weighting function is similar to the kernel function in the SPH method, whereas the kernel function has to satisfy the strict requirement for the SPH type approximation.
367
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
Fig. 1. Particle distribution. (a) uniform distribution, (b) arbitrary distribution. The green circle is the influence domain. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
uj ui ¼ rui ðrj ri Þ u ui j ¼ rui nij
ð5Þ
rj ri
where nij ¼ have
r r
j i . Multiplying by nij on both sides of Eq. (5), we krj ri k
rui ¼
uj ui nij
¼ ðnij nij Þruij
r j ri
ð6Þ
We postulate that Eq. (6) is the relationship between a point r i and one of its neighboring point r j located in the domain Xri . 2 The gradient of u for each particle is unique and will not satisfy Eq. (6) as it is influenced by all the particles within the neighboring domain, the ruij is the micro gradient component generated from neighboring particle j. Therefore, by introducing the weight function on both sides of the Eq. (6) and normalizing the equation, the following formulations can be obtained as
P
ðu u Þ
j i n W krj ri k ij ij P ¼ j2Xr W ij
j2Xri
hP
k2Xrk njk
P
rui :¼ 4
X
k2Xrj
2P ¼4
njk W jk
i
k2Xrj W jk
i
2
rewritten as Eq. (9). However, the components of the matrix A highly depends upon the distance between the center particle and its neighboring particles in the case with the non-uniform distribution as shown in Fig. 1(b). On the other hand, the numerical accuracy would be decreased when the A is replaced by a constant matrix in non-uniform distribution cases.
rui
ð7Þ
3 X ðuj ui Þ 5 4 5 njk njk W jk r j r i nij W ij : j2Xri
P
njk W jk
k2Xrj W jk
31 5
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
2 3 1 4 X ðuj ui Þ nij W ij 5 n0 j2X rj r i
ð8Þ
ri
A
P where n0 ¼ j W 0ij , and n0 has the same definition as it in the traditional MPS method. It is worth noting that for a regular particle within uniformly distribution system as shown in Fig. 1(a), the "P #1 weighted average matrix A ¼
k2Xr
Pj
njk njk W jk
k2Xr
W jk
ð9Þ
ri
Laplacian Operator. To develop the Laplacian operator in terms of the proposed MPS method, the Gauss’s Theorem is exploited with regard to the micro gradient relationship. Recall Gauss’s Theorem with respect to the physical quantity w3 in the domain Xri
Z
V Xr
r wdV ¼ i
Z
@V
w ndS Xr
ð10Þ
i
where n is the normal direction from ri to the boundary @ XX . The same idea of the proposed gradient is utilized in the derivation of r rui by introducing the micro Laplacian r ruij generated from the neighboring point r j , which can be obtained in sense of the Gauss’s Theorem as
r ruij V j ¼ ruij nij Sj
31 2
k2Xrj njk
Ndim X ðuj ui Þ ðrj ri ÞW ij n0 j2X r j r i 2
where ruij corresponds with w in Eq. (10). V j and Sj denote the volume and area of the j in the sense of Gauss’s Theorem, respectively; for instance, a 2-D case is shown in Fig. (2) and the V j and Sj represent the area of the sector and its arc length in the 2-D case, respectively. Note that the dash circle in Fig. (2) is not the influence domain. The ratio c ¼
onal matrix, A ¼ Ndim I; I is identity matrix and N dim equals to dimensional number, which recovers the gradient calculation in the Moving Particle Semi-implicit method, and the Eq. (8) can be 2 The Xri has the same definition as the cut-off domain or influence domain in the standard MPS method.
Vj Sj
in different
problems can be obtained as following: (1) 1-D problem
c¼
in Eq. (8) is a diag-
j
ð11Þ
Vj ¼ r j r i Sj
ð12Þ
(2) 2-D problem
c¼
Vj 1 ¼ r j r i Sj 2
ð13Þ
3 In the Gauss’s Theorem, the physical quantity is a vector or high order tensor, we use the bold w to represent it.
368
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
obtained by combing the corresponding heat flux constitutive model and energy equation. As mentioned before, the C- and Fprocesses heat conduction constitutive model has the capability to better describe the mechanisms of heat conduction process including both the slow Cattaneo-type and fast Fourier-type, and its corresponding temperature equation leads to a generalization of the macro scale theory for heat conduction in solids of the Jeffreys’-type model, Cattaneo model, and the Fourier model for heat conduction; in addition it also inherits a new model in contrast to these mentioned models. In this section, the C–F model and different boundary conditions will be given in general, and the specific derivation of the C–F model can be found in [9]. Define the dimensionless heat conduction number F T as
FT ¼
Fig. 2. Volume and area of the particles in Gauss’s Theorem.
(3) 3-D problem In the 3-D problem, the ratio c is the relationship between the volume of spherical-bottomed conical body and its spherical-bottomed area.
c¼
Vj 1 ¼ r j ri Sj 3
ð14Þ
According to Eq. (6) and Eq. (8), and redefining c ¼ N
1 dim
r j r i ,
KF KF þ KC
ð17Þ
where K F and K C are the heat conductivity due to Fourier process and Cattaneo process, respectively. The one-step temperature formulation of the C–F model is given as follows.
" # @2T @T @2T @ @2T @S ¼ K 2 þ sF T K Cs 2 þ Cs þSþs ; @t @x @t @x2 @t @t C
@T @2T ¼ K 2 þ S; @t @x
for F T ¼ 1 ð18Þ
the micro Laplacian operators of a scalar are as follows:
1
r2 uij ¼ ruij nij
c
31 njk njk W jk uj ui nij 5 P nij rj ri 2 k2Xr W jk
2P
k2Xr j
¼ Ndim 4
ð15Þ
j
It is to be noted that the Laplacian of a scalar is consistent with the derivation of scalar gradient operator, and the Laplacian operator is exploited in the heat conduction problem in terms of r2 T. Therefore, the Laplacian operator in terms of the proposed consistent MPS method can be obtained via micro Laplacian operator: .P 31 P njk njk W jk j2Xri uj ui W ij j2Xr W ij 5 .P i P 2 P rk rj W jk k2Xr W jk W jk
2P
r ui :¼ N dim 4 2
k2Xr j
k2Xrj
j
2P :¼
N dim 4 n0 k
31 X k2Xrj njk njk W jk 5 P uj ui W ij k2Xr W jk j2X j
k2Xrj
ð16Þ
ri
Remark 2.1. (1) In the theoretical derivation of the proposed MPS method, the first order accuracy can be strictly preserved since the first order accuracy of Taylor Expansion is applied; whereas the traditional MPS method is a special case with uniformly distributed particles as mentioned in A ¼ N dim I. In particular, the traditional MPS method is less than the first order approximation in space when the discretization is not uniform. (2) Depending on the Gauss’s Theorem, the proposed Laplacian operator is consistent with the Gradient.
where C is the volumetric heat capacity, S is a heat source, s is the relaxation time, K is the total conductivity, which is sum of the Fourier conductivity K F and the Cattaneo conductivity K C , and K ¼ KF þ KC . It is to be noted that C–F model recovers different heat conduction models via introducing the dimensionless number F T , and when F T 2 ð0; 1Þ the C–F model leads to the Jeffreys’-type model, when F T ¼ 0 the C–F model degenerates to the Cattaneo model; and when F T ¼ 1 the C–F model degenerates to the Fourier model. As such, the governing equations of temperature can be obtained as parabolic (F T ¼ 1) differential equation and hyperbolic (F T ¼ 0) differential equation. In addition, the systems described by the C–F model contains both the first order and second order time derivatives. A generalized i-Integration framework unifies GS4-1 (first-order systems) and GS4-2 (second-order systems) frameworks for the implementation of solving first- or second-order systems without resorting to the individual framework. Therefore, the GS4 i-Integration computational framework is exploited in this work to solve the heat conduction problems governed by the C–F heat conduction model via the consistent MPS formulation. Boundary Conditions. Generally, boundary conditions in the heat conduction problems can be imposed as temperature, heat flux and convection types. In view of the PDEs, temperature is the Dirichlet or Type 1 boundary, heat flux is the Neumann or Type 2 boundary,
2.3. C–F heat conduction constitutive model and boundary conditions The theoretical basis of heat conduction and transport relies on the underlying heat flux constitutive models, and the governing equation describing the heat conduction problems has been
for F T < 1
Fig. 3. Different boundary conditions.
369
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
(2) Type 2: Heat flux The heat flux boundary in the C–F model is determined by the following heat flux expression.
qF ¼ F T K
Fig. 4. Boundary layer.
qC ¼ s
and convection belongs to the Robin or Type 3 boundary, as shown in Fig. 3.
K
K
100Δt 200Δt 300Δt 400Δt 500Δt
0.6
Q_v ¼
v alue ðri 2 Cq Þ
R
Cq q
ð22Þ
ndS
ð23Þ
V Cq
Time 100Δt 200Δt 300Δt 400Δt 500Δt
1.0 0.8 0.6 0.4
0.2
0.2
0.0
ð21Þ
Hence, for each particle located at the boundary Cq , the governing equation considers one more heat source term, and the value is
0.4
0.0
@T ¼ constant @r
1.2
Time
0.8
v alue ðri 2 Cq Þ
Suppose there is a boundary layer exiting along the edge of the geometry in the discretized system as shown in Fig. (4), Let the cross sectional area of the boundary is SCq for the particles located on the boundary layer, and the heat generated per unit volume, per unit time can be obtained as
where T CT ðr CT ; t þ DtÞ is the prescribed boundary values, aj is the particle located in the ghost domain, and bj is the particles symmetrical to particle aj with respect to the boundary particles. (b) Inspired by the application of the prescribed boundaries in the Finite Element Method, the similar approach can be duplicated into the particle based method in the heat conduction problem.
1.0
@T @2T K sF T ¼ constant @r @r@t
If F T ¼ 1
ð19Þ
1.2
ð20Þ
dqC dT ð1 F T ÞK dr dt
Therefore, the heat flux boundary in C–F model can be obtained as follows. If F T < 1
(1) Type 1: Temperature To impose the prescribed boundary temperature into the MPS method, there are two approaches given as follows: (a) A layer of ghost region is added along the boundary of the actual material surfaces underlying the prescribed boundaries. The concept is to treat the actual boundaries as a mirror to complete the influence domain of the boundary particles within one cut-off length. Therefore, the temperature of the particles located in the ghost region at the next time step is defined as:
Tðaj ; t þ DtÞ ¼ 2T CT ðrCT ; t þ DtÞ Tðbj ; tÞ;
dT dr
0.0 0.2
0.4
x
0.6
0.8
1.0
0.0
0.2
(a)
0.4
x
0.6
0.8
1.0
(b)
min s s Fig. 5. Temperature distribution of Case I at different time instants. The solid lines in the plots are exact solution. (a) ðq1 min ; qmax ; qmax 1 ; q1 Þ ¼ ð1;1;1Þ; (b) ðq1 1 ; q1 Þ = (0, 1, 0).
Time Analytical MPS 100Δt 100Δt 200Δt 200Δt 300Δt 300Δt 400Δt 400Δt 500Δt 500Δt
2.5 2.0 1.5
Time Analytical MPS 100Δt 100Δt 200Δt 200Δt 300Δt 300Δt 400Δt 400Δt 500Δt 500Δt
1.4 1.2 1.0 0.8 0.6
1.0
0.4 0.5 0.0 0.0
0.2 0.0 0.2
0.4
x 0.6
(a)
0.8
1.0
0.0
0.2
0.4
x
0.6
0.8
(b)
Fig. 6. Temperature distribution of Case II and III at different time instants. (a) Case II, (b) Case III.
1.0
370
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
If F T ¼ 1
@T q n ¼ K n ¼ bðT 1 Tðr i ; tÞÞðr i 2 CC Þ @r
ð26Þ
where b is convective heat transfer coefficient and T 1 is the temperature of the surrounding medium. Hence, for each particle located at the boundary CC , the governing equation considers one more term, and the value is
Sq ðr i ; tÞ ¼
(a)
(b)
q nSCq 1 ¼ qðri ; tÞ ni V Cq D
ð24Þ
where ni is the normal direction of the boundary and D is the thickness of the boundary layer, D ¼ 2r and r is the radius of particles. (3) Type 3: Convection The convection is the heat transfer between the surface of the body and the surrounding medium. By calculating the heat transfer in a boundary layer as Eq. (24), the convection boundaries can be imposed in the material. The heat transfer due to the convection with regards to each particle located at the boundaries CC is given as follows: If F T < 1
! @T @2T n K sF T q n ¼ K @r @r@t _ i ; tÞ ðri 2 CC Þ ¼ bðT 1 Tðri ; tÞÞ þ bs T_ 1 Tðr
ð27Þ
3. Numerical results and discussion
Fig. 7. Schematic of the problem. (a) Annulus geometry (b) particle discretized system.
Sq ðri ; tÞ ¼
1 qðr i ; tÞ ni D
ð25Þ
1.0
In this section, the proposed consistent MPS method is exploited to simulate several heat conduction problems with different models and boundary effects under the framework of the C–F model and the unified GS4 time integration framework. The specific objectives are as follows: (1) Analyze the numerical simulation of different heat conduction models by the proposed MPS method. A simple 1-D heat conduction problem is exploited to test the proposed consistent MPS method as well as the boundary implementations mentioned in the preceding sections. Under the framework of C–F model, the numerical simulations of different heat conduction models are investigated. The numerical results of the proposed MPS method are compared with the exact solutions. (2) Investigate the ability of the proposed MPS method. A 2-D heat conduction problem with annulus geometry is utilized and different boundary conditions are imposed. The numerical results show the proposed MPS method has the ability to capture the physics.
1.0 Temperature 0 12.50 25.00 37.50 50.00 62.50 75.00 87.50 100.0
0.8 0.6 0.4 0.2 0.0
dT/dt 0.8
0 3338 6675 1.001E4 1.335E4 1.669E4 2.003E4 2.336E4 2.670E4
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
(a)
0.6
0.8
1.0
(b)
1.0
Temperature 0 12.50 25.00 37.50 50.00 62.50 75.00 87.50 100.0
0.8 0.6 0.4 0.2 0.0
1.0 dT/dt 0.8
0 1188 2375 3563 4750 5938 7125 8313 9500
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
(c)
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(d)
Fig. 8. The comparison of different time steps. (a) Temperature distribution at t ¼ 200Dt, (b) temperature changing rate distribution at t ¼ 200Dt, (c) temperature distribution at t ¼ 600Dt, (d) temperature changing rate distribution at t ¼ 600Dt.
371
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
3.1. 1-D heat conduction with different boundary effects In this section, the proposed MPS method is validated and verified by the heat conduction in a 1-D slab, and three different heat conductivity numbers F T are introduced to represent different governing equations, including the first order parabolic (F T ¼ 1) equation and second order hyperbolic (F T ¼ 0) equation, and different types of boundary conditions are imposed in these numerical examples as well. To simplify the mathematical problem, the C and s are set as 1, and the heat source S is not involved in this problem, and the simplified C–F model is given as follows.
" # @ 2 T @T @2T @ @2T ; ¼ K 2 þ FT K þ @x @t @x2 @t2 @t @T @2T ¼K 2; @t @x
for F T < 1 ð28Þ
for F T ¼ 1
To validate and verify the proposed MPS method, the numerical examples used in this part are described as follows: The dimension of this problem is 0 6 x 6 1 and 0 6 t 6 1, and initial condition Tðx; 0Þ ¼ 0. The slab is discretized by 500 material particles and the influence domain has the cut-off radius R ¼ 1:05r, where r is the radius of each particle. (1) Case I: F T =0, Cattaneo Model The left boundary is imposed as Type 1: constant temperature Tð0; tÞ ¼ 1 during the entire simulation duration, and the Type 2 boundary: constant flux K @T j ¼ 1 is imposed @x x¼1 on the right. As shown in Fig. 5, both analytical and the improved MPS results of the temperature distribution are presented at different time instants; and they are in close agreement. Cattaneo–Vernotte equation (Telegraph equation) is degenerated from the C–F model and the numerical results illustrate the wavelike heat conduction. The unphys1.0
Temperature 0.5000 12.94 25.38 37.81 50.25 62.69 75.13 87.56 100.0
0.8 0.6 0.4 0.2 0.0
ical oscillatory behavior on the left side of the discontinuity max s in Fig. 5(a) with GS4 parameters (qmin 1 ; q1 ; q1 )=(1, 1, 1), is decreased by controlling the numerical dissipation s ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) due to the ability of the unified GS4 integration family. Comparing with Fig. 5(a), Fig. 5(b) shows the more physical realistic results. It is worth noting that the results presented in Fig. 5, is purely first order accuracy without any special treatments involving in the proposed approach. Usually, in the traditional methods, such as FEM, FVM or FD, the manners to deal with the discontinuity are through introducing the artificial viscosity, using upwind schemes or considering the high order approximation. These theories and techniques in the proposed MPS method and the other particle-based method would be a future point for our research. (2) Case II: F T =1, Fourier model The left boundary is imposed as Type 2: constant flux @T j ¼ 1 during the entire simulation duration, and the @x x¼0 Type 3 boundary:convection boundary is applied on the right by imposing the environment temperature T 12 jx¼1 ¼ 1 and b ¼ 1. The numerical simulation is carried s on with GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) and both the analytical and numerical results are shown in Fig. 6(a). (3) Case III: F T ¼ 0:1, Jeffreys Model The left boundary is imposed as Type 2: constant flux @T j ¼ 1, and the Type 3 boundary: convection boundary @x x¼0 is applied on the right by imposing the temperature of the surface surrounding the body T 12 jx¼1 ¼ 1 and b ¼ 1. For this s case, GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) is used in the simulation; the influence of constant heat flux and convective type heat flux can be observed in Fig. 6(b), and the numerical results are in good agreement with the analytical results.
1.0 dT/dt 0 3338 6675 1.001E4 1.335E4 1.669E4 2.003E4 2.336E4 2.670E4
0.8 0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
(a)
0.6
0.8
1.0
(b)
1.0
Temperature
0.8
9.500 20.81 32.13 43.44 54.75 66.06 77.38 88.69 100.0
0.6 0.4 0.2 0.0
1.0 dT/dt 0.8
0 1269 2538 3806 5075 6344 7613 8881 1.015E4
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
(c)
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(d)
Fig. 9. The comparison of different time steps. (a) Temperature distribution at t ¼ 200Dt, (b) temperature changing rate distribution at t ¼ 200Dt, (c) temperature distribution at t ¼ 600Dt, (d) temperature changing rate distribution at t ¼ 600Dt.
372
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
It is worth noting that Fig. 6 also presents that the different heat conduction mechanism leads to different physics with the same boundary effects, and due to the feasibility of the C–F model, these physics can be captured easily. 3.2. 2-D heat conduction of annulus structure A 2-D annulus of isotropic material under thermal loads with different boundary conditions, shown in Fig. 7(a), is exploited to investigate the proposed MPS approximation comprehensively, and its discretized system is shown in Fig. 7(b) as well. The inner and outer circles’ radius of the annulus are R1 ¼ 0:15 m and R2 ¼ 0:5 m, respectively, and the thickness of the annulus is 1 m. Its specific heat capacity, thermal conductivity and mass density are specified as cv ¼ 502:4 J/kg-K, k ¼ 16:3 W/m K and q ¼ 7930 kg/m3, respectively. The heat conduction model in this part is focused on the Fourier heat conduction by selecting the heat conductivity number F T ¼ 1 in the C–F model. The following problems are presented as numerical examples. (1) Case I: Two Essential Boundary Conditions Initial condition: Tðx; 0Þ ¼ 0 Boundary condition: T ðx; t ÞjCInner ¼ 0; T ðx; t ÞjCOuter ¼ 100 (2) Case II: One Essential and One Natural Boundary Condition Initial condition: Tðx; 0Þ ¼ 0 Boundary condition: @T ðx;t Þ jCInner @x
¼ 1000; T ðx; tÞjCOuter ¼ 100
Case I The entire domain is discretized by 1220 material parti6
cles and the time step size, Dt ¼ 5 10 s. The influence domain has the cut-off radius R ¼ 1:2r, where r is the radius of one particle; s and GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) is used in the sim_ ulation. The temperature distribution (T) and its changing rate (T) distribution at different time instants are illustrated in Fig. 8. The essential boundary conditions are imposed on the outer and inner edges of the annulus, therefore the temperature changing rate maintains as 0 as the time goes, which is shown in Fig. 8(b) and (d). As a result of the hot boundaries, initial cold temperatures in the system domain increases gradually as the time evolves. However, the temperature changing speed decreases as the temperature distribution within the whole domain turns smooth. Fig. 8 also demonstrates the ability of the proposed MPS in simulating the system with special geometry, and the results preserve similar shapes as that of the system boundary. Case II In this case, the natural boundary is imposed in the annulus’ inner circle. The time integration, spatial discretization and the radius of the influence domain are the same as those of the Case I. _ distribuThe temperature distribution (T) and its changing rate (T) tion at different time instants are illustrated in Fig. 9. Because of the heat flux energy is applied in the inner boundary, the temperature changing speed is time dependent, and the temperature of the entire system is increased by both the hot boundary and heat flux energy.
4. Conclusions A consistent MPS-based particle method has been presented to solve the first- and second-order transient hyperbolic and parabolic systems arising in the study of heat conduction in solids. The proposed MPS-based particle method has been proved to preserve the consistency between the Gradient approximation and the Laplacian approximation; this is in contrast to the relationship between them which does not satisfy the Gauss’s Theorem in the conventional MPS method. Besides, the first order accuracy in space is maintained by the matrix A in the proposed MPS method in the case with arbitrary and non-uniform particle discretization. The approach of implementing different kinds of boundary conditions (Type 1 Dirichlet Boundary condition, Type 2 Neumann boundary condition and Type 3 Robin boundary condition) was validated and verified as well. A theoretical model of solid state heat transfer, the C–F model, was used for illustrating various heat conduction mechanisms, including diffusive, wavelike and ballistic transport, and a unified time integration, GS-4 i-Integration, is exploited to solve the C–F model without the limitation of the time order of governing equations in time. Acknowledgments Acknowledgment is due to the Minnesota Supercomputing Institute (MSI) for computer grants. Tao Xue acknowledges the support by the China Scholarship Council and Nanjing University of Science and Technology. References [1] S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Eng. 123 (3) (1996) 421–434. [2] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [3] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375–389. [4] S. Koshizuka, A. Nobe, Y. Oka, Numerical analysis of breaking waves using the moving particle semi-implicit method, Int. J. Numer. Meth. Fluids 26 (7) (1998) 751–769. [5] A. Khayyer, H. Gotoh, Enhancement of stability and accuracy of the moving particle semi-implicit method, J. Comput. Phys. 230 (8) (2011) 3093–3118. [6] J.B.J. baron Fourier, The Analytical Theory of Heat, The University Press, 1878. [7] M. Carlo Cattaneo, Sulla conduzione de calore, Atti Semin Mat Fis Della Univ Modena 3(3). [8] D.D. Joseph, L. Preziosi, Heat waves, Rev. Mod. Phys. 61 (1) (1989) 41. [9] X. Zhou, K.K. Tamma, C.V. Anderson, On a new C-and F-processes heat conduction constitutive model and the associated generalized theory of dynamic thermoelasticity, J. Therm. Stresses 24 (6) (2001) 531–564. [10] V. Wheeler, S. Masuri, K. Tamma, X. Zhou, M. Sellier, On the applicability of an isochronous integration framework for parabolic/hyperbolic heat conduction type problems, Numer. Heat Transfer, Part A: Appl. 62 (5) (2012) 372–392. [11] M. Shimada, K.K. Tamma, Explicit time integrators and designs for first/second order linear transient systems, Encycl. Therm. Stresses 3 (2013) 1524–1530. [12] M. Shimada, K.K. Tamma, Implicit time integrators and designs for first-/ second-order linear transient systems, Encycl. Therm. Stresses 5 (2013) 2387– 2396. [13] M. Shimada, S. Masuri, K.K. Tamma, Isochronous explicit time integration framework: illustration to thermal stress problems involving both first-and second-order transient systems, J. Therm. Stresses 37 (9) (2014) 1066–1079.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A consistent Moving Particle System Simulation method: Applications to parabolic/hyperbolic heat conduction type problems Tao Xue a,b, Kumar K. Tamma b,⇑, Xiaobing Zhang a a b
School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing, Jiangsu 210094, PR China Department of Mechanical Engineering, University of Minnesota-Twin Cities, 111 Church St. SE, Minneapolis, MN 55455, United States
a r t i c l e
i n f o
Article history: Received 8 March 2016 Received in revised form 5 May 2016 Accepted 6 May 2016 Available online 1 June 2016 Keywords: MPS Particle based method Parabolic/hyperbolic heat conduction GS4 i-Integration
a b s t r a c t In this paper, a consistent MPS (Moving Particle System Semi-implicit/Moving Particle System Simulation) method is established emanating from Taylor Expansion and Gauss’s Theorem. The proposed MPS method preserves the consistency between the Gradient and Laplacian approximations, and also has the ability to recover the conventional MPS method which lacks consistency. In addition, the proposed method preserves the first order accuracy in space whereas the traditional MPS method loses the accuracy in the case with arbitrary and non-uniform particle discretization. The applicability of the proposed method for the analysis of first- and second-order transient hyperbolic and parabolic systems is illustrated via a recently developed C–F model (Cattaneo–Fourier model) of heat conduction that has been proven to govern thermal transport in different scales, and also capture the various heat conduction physical phenomena. The approaches of imposing different boundary conditions in the heat conduction problems are investigated comprehensively as well. In addition, a unified GS4 i-Integration framework is exploited in the numerical simulations involving first- and second-order systems in time. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The MPS (Moving Particle System Semi-implicit/Moving Particle System Simulation) method is a particle based method under the description of Lagrangian meshless method originally proposed for simulations of incompressible fluid dynamics with the objective to capture the free surface [1]. The convenience of MPS method as well as the other particle based meshless methods such as Smooth Particle Hydrodynamics (SPH) method [2,3], is to provide approximations to the strong form of the partial differential equations (PDEs) via using Gradient operator and Laplacian operator. Nevertheless, the (original) MPS method offers the gradient of a scalar value which is mainly used for calculation of the pressure gradient in the Navier–Stokes equation. On the other hand, the Laplacian operator is designed with regards to a vector value which is fluid velocity; consequently this is not consistent with the gradient operator from the perspective of Gauss’s Theorem and basic differential calculus principles. During the past several decades, the MPS method has been widely implemented in various engineering fields. Koshizuka [1,4] established the basic formulations of the MPS method and applied ⇑ Corresponding author. E-mail addresses: [email protected] (T. (K.K. Tamma), [email protected] (X. Zhang).
Xue),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.05.020 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
[email protected]
it to dam break problems for better capturing the free surface. The enhanced MPS method has been developed to improve the stability and accuracy of the MPS methods [5]. One major challenge of the MPS-type simulation in solving PDEs is the boundary condition implementations, which also exists in the other particle based methods. The issue of boundaries in the MPS applications of fluid dynamics is resolved by adding several layers of wall particles, and these wall particles are also used to prevent the fluid particles from penetrating into the wall. However, the approach of imposing boundary conditions, including Dirichlet boundary conditions and Neumann boundary conditions, has not been fully investigated to-date in the MPS method from the perspective of solving a general PDE. Heat Conduction The theoretical basis of heat conduction and transport relies heavily on the heat flux constitutive models and the fusion of the mechanisms of heat transport with regard to heat conduction. Various heat flux models have been proposed for decades and different heat flux models lead to different types of differential equation governing the heat conduction phenomenon in the solid state [6–8]. To better govern the physics of the heat conduction processes and predict the material properties, a novel C-(Cattaneo-type) and F-(Fourier-type) processes conduction constitutive [9,10] has been developed from the physics of the Boltzmann Transport Equation; and the model acknowledges the notion of the simultaneous coexistence of both slow Cattaneo-type C-processes and fast Fourier-type F-processes in the process of
366
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
heat conduction. Different types of heat flux models, lead to the different heat conduction governing equations and represent different mathematical and physical characteristics. In the C–F model, different heat conduction model numbers F T yield different types of PDEs, including the hyperbolic (F T ¼ 0) and parabolic (F T ¼ 1) heat conduction type problems. Hence, the heat conduction C–F model provides us with a comprehensive mathematical example to test and verify the MPS formulation. The Gradient and Laplacian operators can be directly approximated by the MPS method with particle discretization, and the particles carry the physical entity of temperature and transfer it to the neighboring particles. There have been only a limited or scant literature for solving the heat conduction and transport problems via using the MPS method, and the problems involving different types of boundary effects have not been thoroughly investigated. Therefore, the C–F heat conduction model with different type of boundary conditions offers a numerical environment for improving and testing the MPS method, especially from the aspect of the implementation of the boundaries in the PDEs approximated by the MPS method. Time Integration The particle based methods require neighboring particle detecting process for each particle and it turns out to be the very expensive in practical simulations, especially for large scale systems involving large number of particles. Under this general scenario of computational efficiency, the use of explicit time integration schemes with a minimum of second order accuracy in time are recommended in contrast to the first order explicit Euler-based methods. A unified time integration framework, under the umbrella of the so-called generalized single-step single-solve (GSSSS) family of algorithms, including not only implicit, but also explicit algorithms (the representations can be also written in a linear three-step form), are proposed in [11,12]. This simulation framework has been originally developed in the field of computational structural dynamics. It encompasses most existing second order time accurate algorithms, such as Newmark, Central Difference method and HH-a, and more novel and optimal algorithms and advances. Additionally the implementation to first order systems such as heat conduction problem has been developed and proved to retain the same advantages via an isochronous integration [i-Integration] [13] process within the GSSSS framework. In order to improve and extend the traditional MPS method, the main objectives of this paper focus upon the following: (1) To enhance the MPS method to heat conduction problems, an improved consistent MPS-based method is presented depending upon the established MPS Gradient and Laplacian operators. Via the C–F heat flux constitutive model, the consistent MPS-based method is tested and implemented to illustrate the various heat conduction mechanisms in both the parabolic and hyperbolic differential equations. (2) To investigate the implementation of boundaries in the consistent MPS method, as well as other particle-based method, the approach of imposing different types of boundary effects is presented in detail, and these include imposed temperature (Type 1), heat flux (Type 2), convection (Type 3) and the like. (3) To improve the computational efficiency and the numerical compatibility in terms of the developed MPS based method, the second order time accurate unified GS-4 i-Integration framework is implemented within the numerical simulations with regard to both the first- and second- order transient systems.
of the traditional MPS method. Specifically, the terms related to the Gradient and Laplacian of temperature in heat conduction will be developed from Taylor Series, and the consistency between these two operators is strictly preserved via using Gauss’s Theorem; this is in contrast to past practices. 2.1. Standard MPS method In the MPS method the gradient and the Laplacian operators are formulated as:
hr/ii ¼ D
E
r2 /
i
Ndim X /j /i rj ri W ij n0 j–i jjrj r i jj2 ¼
2Ndim X ð/ /i ÞW ij n0 k0 j–i j
ð1Þ
where / is a physical quantity, Ndim is dimensional number of the simulating space, r is coordinate vector, the subscripts i and j are the targeted particle and its neighboring particles, W ij is the weighting function,1 and n0 denotes the particle number density. The formulations of W ij ; n0 and k0 are listed as follows:
W ij ¼
R=rij 1; 0 6 r ij < R
ð2Þ
R6r
0;
where r ij ¼ r j r i ; R is the radius of the influence domain, and
n0 ¼
X j–i
P W 0ij ;
k0 ¼
0 0 j–i W ij r j
P
2 r0i
0 j–i W ij
ð3Þ
It is worth noting that the physical quantity / is a scalar value in the Gradient operator and a vector value in the Laplacian operator; specifically, the objectives of these two operators are to formulate pressure and velocity in the Navier–Stokes equation which the original MPS was developed for. Therefore, these operators are not consistent with each other and the relationship does not satisfy the divergence theorem. In addition, n0 and k0 are obtained by the initial weighting function, W 0ij . 2.2. Consistent MPS-based method We herein present a consistent formulation. The MPS method discretizes the domain as a collection of particles in such a way the physical quantities at each of the particles are affected by its neighboring particles. Consider a domain in Euclidean space, E. Let r i 2 X0 E denotes the geometric position vector of a certain position in the domain and r j 2 X0 E denotes another point which is close to position r i . Therefore, any physical quantity at each point of space and time can be defined as uðri ; tÞ : X0 I ! E. To simplify the formulations, the physical quantity uðr i ; tÞ is replaced by ui in the rest of this paper, the subscript i represents particle label and the subscript, ij, represents the relationship between particle i and j. Gradient Operator. In mathematics, a Taylor series is a representation of a function as an infinite sum of terms that are calculated from the values of the function’s derivatives at a single point. As mentioned above, the first order Taylor Series of a general physical quantity uj about r i at a fixed time t yields
uj ¼ ui þ ruij ðrj ri Þ þ Oðkrj ri k2 Þ
ð4Þ
In the view of Eq. (4), omitting the higher order terms, we have
2. Theory This section recapitulates the main idea used for the derivations of the Gradient and Laplacian operators of a field property in terms
1 The weighting function is similar to the kernel function in the SPH method, whereas the kernel function has to satisfy the strict requirement for the SPH type approximation.
367
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
Fig. 1. Particle distribution. (a) uniform distribution, (b) arbitrary distribution. The green circle is the influence domain. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
uj ui ¼ rui ðrj ri Þ u ui j ¼ rui nij
ð5Þ
rj ri
where nij ¼ have
r r
j i . Multiplying by nij on both sides of Eq. (5), we krj ri k
rui ¼
uj ui nij
¼ ðnij nij Þruij
r j ri
ð6Þ
We postulate that Eq. (6) is the relationship between a point r i and one of its neighboring point r j located in the domain Xri . 2 The gradient of u for each particle is unique and will not satisfy Eq. (6) as it is influenced by all the particles within the neighboring domain, the ruij is the micro gradient component generated from neighboring particle j. Therefore, by introducing the weight function on both sides of the Eq. (6) and normalizing the equation, the following formulations can be obtained as
P
ðu u Þ
j i n W krj ri k ij ij P ¼ j2Xr W ij
j2Xri
hP
k2Xrk njk
P
rui :¼ 4
X
k2Xrj
2P ¼4
njk W jk
i
k2Xrj W jk
i
2
rewritten as Eq. (9). However, the components of the matrix A highly depends upon the distance between the center particle and its neighboring particles in the case with the non-uniform distribution as shown in Fig. 1(b). On the other hand, the numerical accuracy would be decreased when the A is replaced by a constant matrix in non-uniform distribution cases.
rui
ð7Þ
3 X ðuj ui Þ 5 4 5 njk njk W jk r j r i nij W ij : j2Xri
P
njk W jk
k2Xrj W jk
31 5
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
2 3 1 4 X ðuj ui Þ nij W ij 5 n0 j2X rj r i
ð8Þ
ri
A
P where n0 ¼ j W 0ij , and n0 has the same definition as it in the traditional MPS method. It is worth noting that for a regular particle within uniformly distribution system as shown in Fig. 1(a), the "P #1 weighted average matrix A ¼
k2Xr
Pj
njk njk W jk
k2Xr
W jk
ð9Þ
ri
Laplacian Operator. To develop the Laplacian operator in terms of the proposed MPS method, the Gauss’s Theorem is exploited with regard to the micro gradient relationship. Recall Gauss’s Theorem with respect to the physical quantity w3 in the domain Xri
Z
V Xr
r wdV ¼ i
Z
@V
w ndS Xr
ð10Þ
i
where n is the normal direction from ri to the boundary @ XX . The same idea of the proposed gradient is utilized in the derivation of r rui by introducing the micro Laplacian r ruij generated from the neighboring point r j , which can be obtained in sense of the Gauss’s Theorem as
r ruij V j ¼ ruij nij Sj
31 2
k2Xrj njk
Ndim X ðuj ui Þ ðrj ri ÞW ij n0 j2X r j r i 2
where ruij corresponds with w in Eq. (10). V j and Sj denote the volume and area of the j in the sense of Gauss’s Theorem, respectively; for instance, a 2-D case is shown in Fig. (2) and the V j and Sj represent the area of the sector and its arc length in the 2-D case, respectively. Note that the dash circle in Fig. (2) is not the influence domain. The ratio c ¼
onal matrix, A ¼ Ndim I; I is identity matrix and N dim equals to dimensional number, which recovers the gradient calculation in the Moving Particle Semi-implicit method, and the Eq. (8) can be 2 The Xri has the same definition as the cut-off domain or influence domain in the standard MPS method.
Vj Sj
in different
problems can be obtained as following: (1) 1-D problem
c¼
in Eq. (8) is a diag-
j
ð11Þ
Vj ¼ r j r i Sj
ð12Þ
(2) 2-D problem
c¼
Vj 1 ¼ r j r i Sj 2
ð13Þ
3 In the Gauss’s Theorem, the physical quantity is a vector or high order tensor, we use the bold w to represent it.
368
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
obtained by combing the corresponding heat flux constitutive model and energy equation. As mentioned before, the C- and Fprocesses heat conduction constitutive model has the capability to better describe the mechanisms of heat conduction process including both the slow Cattaneo-type and fast Fourier-type, and its corresponding temperature equation leads to a generalization of the macro scale theory for heat conduction in solids of the Jeffreys’-type model, Cattaneo model, and the Fourier model for heat conduction; in addition it also inherits a new model in contrast to these mentioned models. In this section, the C–F model and different boundary conditions will be given in general, and the specific derivation of the C–F model can be found in [9]. Define the dimensionless heat conduction number F T as
FT ¼
Fig. 2. Volume and area of the particles in Gauss’s Theorem.
(3) 3-D problem In the 3-D problem, the ratio c is the relationship between the volume of spherical-bottomed conical body and its spherical-bottomed area.
c¼
Vj 1 ¼ r j ri Sj 3
ð14Þ
According to Eq. (6) and Eq. (8), and redefining c ¼ N
1 dim
r j r i ,
KF KF þ KC
ð17Þ
where K F and K C are the heat conductivity due to Fourier process and Cattaneo process, respectively. The one-step temperature formulation of the C–F model is given as follows.
" # @2T @T @2T @ @2T @S ¼ K 2 þ sF T K Cs 2 þ Cs þSþs ; @t @x @t @x2 @t @t C
@T @2T ¼ K 2 þ S; @t @x
for F T ¼ 1 ð18Þ
the micro Laplacian operators of a scalar are as follows:
1
r2 uij ¼ ruij nij
c
31 njk njk W jk uj ui nij 5 P nij rj ri 2 k2Xr W jk
2P
k2Xr j
¼ Ndim 4
ð15Þ
j
It is to be noted that the Laplacian of a scalar is consistent with the derivation of scalar gradient operator, and the Laplacian operator is exploited in the heat conduction problem in terms of r2 T. Therefore, the Laplacian operator in terms of the proposed consistent MPS method can be obtained via micro Laplacian operator: .P 31 P njk njk W jk j2Xri uj ui W ij j2Xr W ij 5 .P i P 2 P rk rj W jk k2Xr W jk W jk
2P
r ui :¼ N dim 4 2
k2Xr j
k2Xrj
j
2P :¼
N dim 4 n0 k
31 X k2Xrj njk njk W jk 5 P uj ui W ij k2Xr W jk j2X j
k2Xrj
ð16Þ
ri
Remark 2.1. (1) In the theoretical derivation of the proposed MPS method, the first order accuracy can be strictly preserved since the first order accuracy of Taylor Expansion is applied; whereas the traditional MPS method is a special case with uniformly distributed particles as mentioned in A ¼ N dim I. In particular, the traditional MPS method is less than the first order approximation in space when the discretization is not uniform. (2) Depending on the Gauss’s Theorem, the proposed Laplacian operator is consistent with the Gradient.
where C is the volumetric heat capacity, S is a heat source, s is the relaxation time, K is the total conductivity, which is sum of the Fourier conductivity K F and the Cattaneo conductivity K C , and K ¼ KF þ KC . It is to be noted that C–F model recovers different heat conduction models via introducing the dimensionless number F T , and when F T 2 ð0; 1Þ the C–F model leads to the Jeffreys’-type model, when F T ¼ 0 the C–F model degenerates to the Cattaneo model; and when F T ¼ 1 the C–F model degenerates to the Fourier model. As such, the governing equations of temperature can be obtained as parabolic (F T ¼ 1) differential equation and hyperbolic (F T ¼ 0) differential equation. In addition, the systems described by the C–F model contains both the first order and second order time derivatives. A generalized i-Integration framework unifies GS4-1 (first-order systems) and GS4-2 (second-order systems) frameworks for the implementation of solving first- or second-order systems without resorting to the individual framework. Therefore, the GS4 i-Integration computational framework is exploited in this work to solve the heat conduction problems governed by the C–F heat conduction model via the consistent MPS formulation. Boundary Conditions. Generally, boundary conditions in the heat conduction problems can be imposed as temperature, heat flux and convection types. In view of the PDEs, temperature is the Dirichlet or Type 1 boundary, heat flux is the Neumann or Type 2 boundary,
2.3. C–F heat conduction constitutive model and boundary conditions The theoretical basis of heat conduction and transport relies on the underlying heat flux constitutive models, and the governing equation describing the heat conduction problems has been
for F T < 1
Fig. 3. Different boundary conditions.
369
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
(2) Type 2: Heat flux The heat flux boundary in the C–F model is determined by the following heat flux expression.
qF ¼ F T K
Fig. 4. Boundary layer.
qC ¼ s
and convection belongs to the Robin or Type 3 boundary, as shown in Fig. 3.
K
K
100Δt 200Δt 300Δt 400Δt 500Δt
0.6
Q_v ¼
v alue ðri 2 Cq Þ
R
Cq q
ð22Þ
ndS
ð23Þ
V Cq
Time 100Δt 200Δt 300Δt 400Δt 500Δt
1.0 0.8 0.6 0.4
0.2
0.2
0.0
ð21Þ
Hence, for each particle located at the boundary Cq , the governing equation considers one more heat source term, and the value is
0.4
0.0
@T ¼ constant @r
1.2
Time
0.8
v alue ðri 2 Cq Þ
Suppose there is a boundary layer exiting along the edge of the geometry in the discretized system as shown in Fig. (4), Let the cross sectional area of the boundary is SCq for the particles located on the boundary layer, and the heat generated per unit volume, per unit time can be obtained as
where T CT ðr CT ; t þ DtÞ is the prescribed boundary values, aj is the particle located in the ghost domain, and bj is the particles symmetrical to particle aj with respect to the boundary particles. (b) Inspired by the application of the prescribed boundaries in the Finite Element Method, the similar approach can be duplicated into the particle based method in the heat conduction problem.
1.0
@T @2T K sF T ¼ constant @r @r@t
If F T ¼ 1
ð19Þ
1.2
ð20Þ
dqC dT ð1 F T ÞK dr dt
Therefore, the heat flux boundary in C–F model can be obtained as follows. If F T < 1
(1) Type 1: Temperature To impose the prescribed boundary temperature into the MPS method, there are two approaches given as follows: (a) A layer of ghost region is added along the boundary of the actual material surfaces underlying the prescribed boundaries. The concept is to treat the actual boundaries as a mirror to complete the influence domain of the boundary particles within one cut-off length. Therefore, the temperature of the particles located in the ghost region at the next time step is defined as:
Tðaj ; t þ DtÞ ¼ 2T CT ðrCT ; t þ DtÞ Tðbj ; tÞ;
dT dr
0.0 0.2
0.4
x
0.6
0.8
1.0
0.0
0.2
(a)
0.4
x
0.6
0.8
1.0
(b)
min s s Fig. 5. Temperature distribution of Case I at different time instants. The solid lines in the plots are exact solution. (a) ðq1 min ; qmax ; qmax 1 ; q1 Þ ¼ ð1;1;1Þ; (b) ðq1 1 ; q1 Þ = (0, 1, 0).
Time Analytical MPS 100Δt 100Δt 200Δt 200Δt 300Δt 300Δt 400Δt 400Δt 500Δt 500Δt
2.5 2.0 1.5
Time Analytical MPS 100Δt 100Δt 200Δt 200Δt 300Δt 300Δt 400Δt 400Δt 500Δt 500Δt
1.4 1.2 1.0 0.8 0.6
1.0
0.4 0.5 0.0 0.0
0.2 0.0 0.2
0.4
x 0.6
(a)
0.8
1.0
0.0
0.2
0.4
x
0.6
0.8
(b)
Fig. 6. Temperature distribution of Case II and III at different time instants. (a) Case II, (b) Case III.
1.0
370
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
If F T ¼ 1
@T q n ¼ K n ¼ bðT 1 Tðr i ; tÞÞðr i 2 CC Þ @r
ð26Þ
where b is convective heat transfer coefficient and T 1 is the temperature of the surrounding medium. Hence, for each particle located at the boundary CC , the governing equation considers one more term, and the value is
Sq ðr i ; tÞ ¼
(a)
(b)
q nSCq 1 ¼ qðri ; tÞ ni V Cq D
ð24Þ
where ni is the normal direction of the boundary and D is the thickness of the boundary layer, D ¼ 2r and r is the radius of particles. (3) Type 3: Convection The convection is the heat transfer between the surface of the body and the surrounding medium. By calculating the heat transfer in a boundary layer as Eq. (24), the convection boundaries can be imposed in the material. The heat transfer due to the convection with regards to each particle located at the boundaries CC is given as follows: If F T < 1
! @T @2T n K sF T q n ¼ K @r @r@t _ i ; tÞ ðri 2 CC Þ ¼ bðT 1 Tðri ; tÞÞ þ bs T_ 1 Tðr
ð27Þ
3. Numerical results and discussion
Fig. 7. Schematic of the problem. (a) Annulus geometry (b) particle discretized system.
Sq ðri ; tÞ ¼
1 qðr i ; tÞ ni D
ð25Þ
1.0
In this section, the proposed consistent MPS method is exploited to simulate several heat conduction problems with different models and boundary effects under the framework of the C–F model and the unified GS4 time integration framework. The specific objectives are as follows: (1) Analyze the numerical simulation of different heat conduction models by the proposed MPS method. A simple 1-D heat conduction problem is exploited to test the proposed consistent MPS method as well as the boundary implementations mentioned in the preceding sections. Under the framework of C–F model, the numerical simulations of different heat conduction models are investigated. The numerical results of the proposed MPS method are compared with the exact solutions. (2) Investigate the ability of the proposed MPS method. A 2-D heat conduction problem with annulus geometry is utilized and different boundary conditions are imposed. The numerical results show the proposed MPS method has the ability to capture the physics.
1.0 Temperature 0 12.50 25.00 37.50 50.00 62.50 75.00 87.50 100.0
0.8 0.6 0.4 0.2 0.0
dT/dt 0.8
0 3338 6675 1.001E4 1.335E4 1.669E4 2.003E4 2.336E4 2.670E4
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
(a)
0.6
0.8
1.0
(b)
1.0
Temperature 0 12.50 25.00 37.50 50.00 62.50 75.00 87.50 100.0
0.8 0.6 0.4 0.2 0.0
1.0 dT/dt 0.8
0 1188 2375 3563 4750 5938 7125 8313 9500
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
(c)
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(d)
Fig. 8. The comparison of different time steps. (a) Temperature distribution at t ¼ 200Dt, (b) temperature changing rate distribution at t ¼ 200Dt, (c) temperature distribution at t ¼ 600Dt, (d) temperature changing rate distribution at t ¼ 600Dt.
371
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
3.1. 1-D heat conduction with different boundary effects In this section, the proposed MPS method is validated and verified by the heat conduction in a 1-D slab, and three different heat conductivity numbers F T are introduced to represent different governing equations, including the first order parabolic (F T ¼ 1) equation and second order hyperbolic (F T ¼ 0) equation, and different types of boundary conditions are imposed in these numerical examples as well. To simplify the mathematical problem, the C and s are set as 1, and the heat source S is not involved in this problem, and the simplified C–F model is given as follows.
" # @ 2 T @T @2T @ @2T ; ¼ K 2 þ FT K þ @x @t @x2 @t2 @t @T @2T ¼K 2; @t @x
for F T < 1 ð28Þ
for F T ¼ 1
To validate and verify the proposed MPS method, the numerical examples used in this part are described as follows: The dimension of this problem is 0 6 x 6 1 and 0 6 t 6 1, and initial condition Tðx; 0Þ ¼ 0. The slab is discretized by 500 material particles and the influence domain has the cut-off radius R ¼ 1:05r, where r is the radius of each particle. (1) Case I: F T =0, Cattaneo Model The left boundary is imposed as Type 1: constant temperature Tð0; tÞ ¼ 1 during the entire simulation duration, and the Type 2 boundary: constant flux K @T j ¼ 1 is imposed @x x¼1 on the right. As shown in Fig. 5, both analytical and the improved MPS results of the temperature distribution are presented at different time instants; and they are in close agreement. Cattaneo–Vernotte equation (Telegraph equation) is degenerated from the C–F model and the numerical results illustrate the wavelike heat conduction. The unphys1.0
Temperature 0.5000 12.94 25.38 37.81 50.25 62.69 75.13 87.56 100.0
0.8 0.6 0.4 0.2 0.0
ical oscillatory behavior on the left side of the discontinuity max s in Fig. 5(a) with GS4 parameters (qmin 1 ; q1 ; q1 )=(1, 1, 1), is decreased by controlling the numerical dissipation s ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) due to the ability of the unified GS4 integration family. Comparing with Fig. 5(a), Fig. 5(b) shows the more physical realistic results. It is worth noting that the results presented in Fig. 5, is purely first order accuracy without any special treatments involving in the proposed approach. Usually, in the traditional methods, such as FEM, FVM or FD, the manners to deal with the discontinuity are through introducing the artificial viscosity, using upwind schemes or considering the high order approximation. These theories and techniques in the proposed MPS method and the other particle-based method would be a future point for our research. (2) Case II: F T =1, Fourier model The left boundary is imposed as Type 2: constant flux @T j ¼ 1 during the entire simulation duration, and the @x x¼0 Type 3 boundary:convection boundary is applied on the right by imposing the environment temperature T 12 jx¼1 ¼ 1 and b ¼ 1. The numerical simulation is carried s on with GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) and both the analytical and numerical results are shown in Fig. 6(a). (3) Case III: F T ¼ 0:1, Jeffreys Model The left boundary is imposed as Type 2: constant flux @T j ¼ 1, and the Type 3 boundary: convection boundary @x x¼0 is applied on the right by imposing the temperature of the surface surrounding the body T 12 jx¼1 ¼ 1 and b ¼ 1. For this s case, GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) is used in the simulation; the influence of constant heat flux and convective type heat flux can be observed in Fig. 6(b), and the numerical results are in good agreement with the analytical results.
1.0 dT/dt 0 3338 6675 1.001E4 1.335E4 1.669E4 2.003E4 2.336E4 2.670E4
0.8 0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
(a)
0.6
0.8
1.0
(b)
1.0
Temperature
0.8
9.500 20.81 32.13 43.44 54.75 66.06 77.38 88.69 100.0
0.6 0.4 0.2 0.0
1.0 dT/dt 0.8
0 1269 2538 3806 5075 6344 7613 8881 1.015E4
0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
(c)
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(d)
Fig. 9. The comparison of different time steps. (a) Temperature distribution at t ¼ 200Dt, (b) temperature changing rate distribution at t ¼ 200Dt, (c) temperature distribution at t ¼ 600Dt, (d) temperature changing rate distribution at t ¼ 600Dt.
372
T. Xue et al. / International Journal of Heat and Mass Transfer 101 (2016) 365–372
It is worth noting that Fig. 6 also presents that the different heat conduction mechanism leads to different physics with the same boundary effects, and due to the feasibility of the C–F model, these physics can be captured easily. 3.2. 2-D heat conduction of annulus structure A 2-D annulus of isotropic material under thermal loads with different boundary conditions, shown in Fig. 7(a), is exploited to investigate the proposed MPS approximation comprehensively, and its discretized system is shown in Fig. 7(b) as well. The inner and outer circles’ radius of the annulus are R1 ¼ 0:15 m and R2 ¼ 0:5 m, respectively, and the thickness of the annulus is 1 m. Its specific heat capacity, thermal conductivity and mass density are specified as cv ¼ 502:4 J/kg-K, k ¼ 16:3 W/m K and q ¼ 7930 kg/m3, respectively. The heat conduction model in this part is focused on the Fourier heat conduction by selecting the heat conductivity number F T ¼ 1 in the C–F model. The following problems are presented as numerical examples. (1) Case I: Two Essential Boundary Conditions Initial condition: Tðx; 0Þ ¼ 0 Boundary condition: T ðx; t ÞjCInner ¼ 0; T ðx; t ÞjCOuter ¼ 100 (2) Case II: One Essential and One Natural Boundary Condition Initial condition: Tðx; 0Þ ¼ 0 Boundary condition: @T ðx;t Þ jCInner @x
¼ 1000; T ðx; tÞjCOuter ¼ 100
Case I The entire domain is discretized by 1220 material parti6
cles and the time step size, Dt ¼ 5 10 s. The influence domain has the cut-off radius R ¼ 1:2r, where r is the radius of one particle; s and GS4 parameters ðq1 min ; qmax 1 ; q1 Þ = (0, 1, 0) is used in the sim_ ulation. The temperature distribution (T) and its changing rate (T) distribution at different time instants are illustrated in Fig. 8. The essential boundary conditions are imposed on the outer and inner edges of the annulus, therefore the temperature changing rate maintains as 0 as the time goes, which is shown in Fig. 8(b) and (d). As a result of the hot boundaries, initial cold temperatures in the system domain increases gradually as the time evolves. However, the temperature changing speed decreases as the temperature distribution within the whole domain turns smooth. Fig. 8 also demonstrates the ability of the proposed MPS in simulating the system with special geometry, and the results preserve similar shapes as that of the system boundary. Case II In this case, the natural boundary is imposed in the annulus’ inner circle. The time integration, spatial discretization and the radius of the influence domain are the same as those of the Case I. _ distribuThe temperature distribution (T) and its changing rate (T) tion at different time instants are illustrated in Fig. 9. Because of the heat flux energy is applied in the inner boundary, the temperature changing speed is time dependent, and the temperature of the entire system is increased by both the hot boundary and heat flux energy.
4. Conclusions A consistent MPS-based particle method has been presented to solve the first- and second-order transient hyperbolic and parabolic systems arising in the study of heat conduction in solids. The proposed MPS-based particle method has been proved to preserve the consistency between the Gradient approximation and the Laplacian approximation; this is in contrast to the relationship between them which does not satisfy the Gauss’s Theorem in the conventional MPS method. Besides, the first order accuracy in space is maintained by the matrix A in the proposed MPS method in the case with arbitrary and non-uniform particle discretization. The approach of implementing different kinds of boundary conditions (Type 1 Dirichlet Boundary condition, Type 2 Neumann boundary condition and Type 3 Robin boundary condition) was validated and verified as well. A theoretical model of solid state heat transfer, the C–F model, was used for illustrating various heat conduction mechanisms, including diffusive, wavelike and ballistic transport, and a unified time integration, GS-4 i-Integration, is exploited to solve the C–F model without the limitation of the time order of governing equations in time. Acknowledgments Acknowledgment is due to the Minnesota Supercomputing Institute (MSI) for computer grants. Tao Xue acknowledges the support by the China Scholarship Council and Nanjing University of Science and Technology. References [1] S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Eng. 123 (3) (1996) 421–434. [2] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [3] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375–389. [4] S. Koshizuka, A. Nobe, Y. Oka, Numerical analysis of breaking waves using the moving particle semi-implicit method, Int. J. Numer. Meth. Fluids 26 (7) (1998) 751–769. [5] A. Khayyer, H. Gotoh, Enhancement of stability and accuracy of the moving particle semi-implicit method, J. Comput. Phys. 230 (8) (2011) 3093–3118. [6] J.B.J. baron Fourier, The Analytical Theory of Heat, The University Press, 1878. [7] M. Carlo Cattaneo, Sulla conduzione de calore, Atti Semin Mat Fis Della Univ Modena 3(3). [8] D.D. Joseph, L. Preziosi, Heat waves, Rev. Mod. Phys. 61 (1) (1989) 41. [9] X. Zhou, K.K. Tamma, C.V. Anderson, On a new C-and F-processes heat conduction constitutive model and the associated generalized theory of dynamic thermoelasticity, J. Therm. Stresses 24 (6) (2001) 531–564. [10] V. Wheeler, S. Masuri, K. Tamma, X. Zhou, M. Sellier, On the applicability of an isochronous integration framework for parabolic/hyperbolic heat conduction type problems, Numer. Heat Transfer, Part A: Appl. 62 (5) (2012) 372–392. [11] M. Shimada, K.K. Tamma, Explicit time integrators and designs for first/second order linear transient systems, Encycl. Therm. Stresses 3 (2013) 1524–1530. [12] M. Shimada, K.K. Tamma, Implicit time integrators and designs for first-/ second-order linear transient systems, Encycl. Therm. Stresses 5 (2013) 2387– 2396. [13] M. Shimada, S. Masuri, K.K. Tamma, Isochronous explicit time integration framework: illustration to thermal stress problems involving both first-and second-order transient systems, J. Therm. Stresses 37 (9) (2014) 1066–1079.