- Email: [email protected]
Generalized heat conduction model involving imperfect thermal contact surface: Application of the GSSSS-1 differential-algebraic equation time integration
International Journal of Heat and Mass Transfer 116 (2018) 889–896
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Generalized heat conduction model involving imperfect thermal contact surface: Application of the GSSSS-1 differential-algebraic equation time integration Tao Xue a, Xiaobing Zhang a,⇑, Kumar K. Tamma b 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, Minneapolis, USA
a r t i c l e
i n f o
Article history: Received 6 May 2017 Received in revised form 19 September 2017 Accepted 21 September 2017 Available online 17 October 2017 Keywords: Thermal contact Imperfect interface Differential algebraic equation Generalized heat conduction model
a b s t r a c t A number of contributions have been made to model thermal contact problems involving perfect/imperfect thermal interfaces. The Fourier-, Cattaneo-, and Jeffreys-types of thermal flux modes have been exploited to govern the heat transfer processes of adjacent regions. However, most of the existing studies consider using only one type of these thermal flux models to describe the thermal physics of the entire system such that the mathematical model may fail to precisely describe the physics within the multidomain system, especially, when the mean free paths of different adjacent domains span in a large range with respect to the characteristic dimension of the material. To circumvent this issue, a generalized heat conduction model, termed C- and F-processes heat conduction model, with the consideration of perfect/ imperfect thermal interface is proposed and formulated under the expression of a first-order time dependent differential-algebraic equation in which the interface conditions are treated as algebraic equations. The unified time integration of GSSSS-1 is extended to solve the resulting differential-algebraic system with Index 2 constraints. The numerical results illustrate that the proposed framework has a better capacity of describing the thermal contacts with/without the thermal resistance of Fourier-Fourier, Fourier-Cattaneo, and other general cases. Ó 2017 Published by Elsevier Ltd.
1. Introduction For the coupling of different scales/materials within a solid, of either multilayered nano-scale films, micro devices or macro systems, the interfaces connecting adjacent domains play an important role in the thermal response of the entire system. In the ideal situation, the adjacent domains contact perfectly such that the temperature and the flux perform no jumps across the interface due to the energy conservation law. The interface of the perfect contacting case is classified as Perfect Interface. However, in the general practical applications, the adjacent domains usually do not contact well and there exist gaps between the contacting domains. In such a way, this Imperfect Interface leads a discontinuity of the temperature while the thermal flux retains a continuity across the interface. Specifically, the relationships with respect to the temperature T and the thermal flux q for the imperfect interfaces [1–3] are given by: ⇑ Corresponding author. E-mail addresses: [email protected] (T. Xue), [email protected] (X. Zhang), [email protected] (K.K. Tamma). https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.081 0017-9310/Ó 2017 Published by Elsevier Ltd.
sTt ¼ Rffqgg n sqt n ¼ 0
ð1Þ
where sTt represents TjC1 TjC2 and ffqgg represents qjC1 þ qjC2 . The coefficient R P 0 represents the thermal (Kapitza)
1 2
resistance, when R ¼ 0 the thermal contact is associated with the perfect interface while R – 0 the thermal contact is associated with the imperfect interface. C1 and C2 are the surfaces in contact and n is the normal direction from C2 to C1 . It is worth noting that Eq. (1) represents the lowly-conducting interface, which is classically modeled using Kapitza’s assumption of thermal resistance. A more general formulation regarding the interface conditions has been proposed in [4], both the lowly-conducting and the highlyconducting interfaces are covered in this general imperfect interface theory. However, the general imperfect interface assumes that the interface is a zero-thickness material and has its own heat transfer process such that it may introduce more numerical problems. For simplicity, the lowly-conducting interface is of interest to this article. Numerous analytical, experimental, and numerical models have been done regarding the thermal contact and the treatments of the
890
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
interface condition. In some thermal contact problems, the adjacent regions are governed by the Fourier’ s heat conduction model. The weak-form based numerical approaches such as FEM, XFEM are applied to capture the strong discontinuity caused by the thermal resistance [5,1,3]. In the engineering applications of thermal wave theory involving thermal interfaces, the thermal contact is usually considered in the framework of the hyperbolic heat conduction model. The thermal resistance was treated as a jump boundary conditions in [6,7]. Additionally, remarkable studies regarding the theoretical analysis with respect to the thermal contact with/without resistance in the framework of hyperbolic model have been done in [8,9]. Note that the numerical approaches applied in the hyperbolic-type system should be able to capture the discontinuities in each domain and also works well under the condition of thermal resistance. The thermal contact has also been extended to the field of bio-engineering, such as the thermal processes in the multi-layer bio-tissues. Both the Fourier-type heat conduction and dual-phase-lag heat conduction model1 have been done with respect to this particular aspect [11,12]. Instead of the linear thermal resistance between the temperature and the flux at the interface, the heat flux across the interface of two layers is assumed to proportionally change with the difference between the fourth powers of the temperatures on the two sides of the interface, which is similar to the radiation heat transfer model. The specific expression of R [8] is as follows:
TjC TjC2 R¼ 1 T4j T4j h C1 C2
ð2Þ
is a thermal contact parameter related to the surface where h roughness [13]. The overview described earlier indicates that the existing works consider implementing only one certain thermal flux model to govern the thermal transfer over the entire multi-domain system, which seems to be not precise and cannot be considered as a complete theory when the mean free paths of different adjacent domains span in a large range with respect to the characteristic dimension of the material. For example, [6] utilized the hyperbolic heat conduction model in all the adjacent regions to formulate the thermal contact of the practical electronic device, which are comprised by conducting film and substrate composite. However, the proper thermal flux models for the film and the subtract are necessary to be determined by the relationship between the mean free path (k) and the characteristic dimension of the material (L). Specifically, 1. When k L, the pure ballistic transport dominates the transport process; 2. When k L, the transport process is governed by a diffusiveballistic nature; 3. When k L, the heat transport is treat to be macroscopic and is governed by a purely diffusive nature. In other words, using a single heat conduction model cannot fulfill the requirement of precisely describing the thermal process over the entire system. In order to achieve it, different types of heat conduction models are coupled by the interface, which turns out to be a numerical challenge for the previous studies. The specific reasons include that: (i) different orders of time-dependent systems are needed to describe the thermal contact between the Fourier-type flux and the Cattaneo-/Jeffreys-type flux; 1 The underlying physics of the dual-phase-lag has been proved to be the same as the Jeffreys’ conduction model. More details can be found in [10].
(ii) the flux models of the Cattaneo-type and the Jeffreys-type involve the partial time derivative and are not explicitly given, which implies that the interface condition of thermal flux is not easy to be fulfilled; (iii) the existing numerical treatments with respect to both the perfect and imperfect interfaces also seem to be a numerical challenge in different approaches, specifically, the interface is usually handled as a boundary condition that are embedded in the time stepping of the governing equation. Consequently, in most of the previous studies the governing equations of each domain have to be evaluated separately and the time stepping strategies to the associated numerical simulations are designed in the framework of staggering methodology, which may have the issue of violating the interface conditions. Therefore, increasing application and necessity regarding the treatment of interface, the modeling of the heat conduction within adjacent regions, and the thermal performance at different scales have motivated the objective of this work to develop a general framework of describing the complex thermal contact physics and a decent approach to deal with the perfect/imperfect interface condition. To circumvent these issues, a two-field (both the temperature and the thermal flux are primary variables) generalized heat conduction model, the C- and F-processes heat conduction model, is utilized to depict the heat conduction of the adjacent regions. The C- and F-processes heat conduction model is originally developed in [14] and is based upon the assumption that the Cattaneo-type thermal flux associated with a finite sound speed and the Fourier-type thermal flux associated with an infinite sound speed co-exist in a heat process simultaneously. Therefore, the C- and F-heat conduction model is acknowledged to cover various heat conduction models. Though the hyperbolicand Jeffreys-types of heat conduction are acknowledged as second-order time systems in terms of single-field formulation (only the temperature is the primary variable), the governing equation of the proposed C-F heat conduction is presented in a first-order two-field form with considering both the temperature and the thermal flux as the primary physical fields of interest. In the framework of the two-field C-and F-model, the thermal flux model in each adjacent region can be switched easily and all the governing equations can be expressed in a unified formulation. Different from the past efforts, we treat the interface condition, Eq. (1), as an algebraic equation such that coupling the interface algebraic equation with the generalized two-field C- and F-heat conduction equation formulates a first-order time dependent Differential Algebraic Equation. The principal contributions emanating from such a generalized framework are summarized as follows: 1. A simple computational framework to model the thermal contact between either the Fourier-type and/or Cattaneo-type and/or Jeffreys-type, without having to resort to the individual governing equation. 2. A perfect preservation of the perfect/imperfect thermal interface condition in a machine error (1E-16). Additionally, we also present the formalism of a unified secondorder accurate time integration, GS4 (generalized single step single solve)[15–17], to solve this first order DAEs system in terms of GS4-1 DAE Index-2. Note that the Index 2 is the highest constraint index existing in the first-order time dependent DAE system and naturally arises from the interface condition in this particular thermal contact problem.
891
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
2. F T 2 ð0; 1Þ: Jeffreys-type Model
2. Theory In this section, we first introduce the fundamental theory of the generalized C- and F-processes model and its associated two-field form. Then, the differential algebraic equation is formulated to describe the heat conduction model involving the perfect/imperfect thermal interface. Last, the GS4-1 DAE Index 2 time integration is derived in detail.
The standard C-and F-heat conduction model was developed from the fundamental Boltzmann Transport Equation under the assumption that the Fourier-type and the Cattaneo-type processes coexist simultaneously in the heat conduction process [14,10]. Compared with the Jeffreys-type heat conduction model (DualPhase-Lag heat conduction), which can only reduce to a Fourierlike heat conduction model with relaxation, the C- and F-model not only explains the underlying theory of the slow processes (Cattaneo-type), fast processes (Fourier-type), and those in between; but also naturally leads to a generalization of the heat conduction in solids of the Jeffreys-, Cattaneo-, and Fourier-type models. The standard C-and F-model is given by:
ð3Þ
where s represents the relaxation time; q; qF , and qC are the total heat flux, the fast heat flux (Fourier-type), and the slow heat flux (Cattaneo-type), respectively. The corresponding heat conductiviF ties of q; qF , and qC are k; kF , and kC , respectively. F T ¼ kCKþk is defined F as a macro-scale heat conduction model number to physically illustrate different types of heat conduction processes, specifically, 1. F T ¼ 0 The C- and F-model degenerates to the Cattaneo-type (hyperbolic) heat conduction processes, which has a finite heat propagation speed. 2. F T ¼ 1 The C- and F-model degenerates to the Fourier-type (parabolic) heat conduction processes, which has an infinite heat propagation speed. 3. F T 2 ð0; 1Þ The C- and F-model degenerates to the Jeffreys-type (parabolic) heat conduction processes, the discontinuities existing in the Jeffreys-type model is smoothed by the diffusion (Fouriertype) effects. It’s worth noting that the acknowledged Jeffreystype heat conduction is also considered as a generalized model, however, it’s not able to obtain the Fourier-type heat conduction. The conservation of energy yields the following
1 @T ¼rqþS a @t
ð4Þ
where a ¼ qkc ; q and c represent the material density and specific heat, respectively; S is a heat source within the material. Substituting the C-and F-heat flux model to Eq. (4), the corresponding generalized single-field form with respect to the one temperature theory are illustrated as follows: 1. F T ¼ 0: Cattaneo-type Model
1 @ T 1 @T 1 @S 2 þ r T þ S þ s ¼ @t s c2T @t2 a @t 2
ð5Þ
ð6Þ
3. F T ¼ 1: Fourier-type Model
1 @T
a @t
¼ r2 T þ S
where cT ¼
2.1. Generalized heat conduction C- and F-model
qF ¼ kF T rT @q qC þ s C ¼ ð1 F T ÞkrT @t q ¼ qF þ qC
1 @ 2 T 1 @T @ r2 T 1 @S 2 þ r T þ s F S þ s ¼ þ T @t @t s c2T @t 2 a @t
paffiffi s
ð7Þ
in the Cattaneo- and Jeffreys-type models.
2.2. Differential algebraic equations of thermal contact based on the Cand F-heat conduction model In this work, instead of the traditional single-field form, we propose to directly utilize the two-field form of the C- and F-model in the aspect of the thermal contact problem. Both the temperature and the thermal flux are treated as the primary variables. The interface conditions are naturally comprised by these two primary variables and considered to be the constraints with respect to the governing equation. To formulate the DAE with Index 2 constraint, the interface conditions are treated as constraint and a Lagrangian Multiplier (k) is introduced to couple the constraints with the generalized C- and F-heat conduction model. The specific formulation is as follows:
M q_ þ Kq þ BT k ¼ F UðqÞ ¼ 0
ð8Þ
h iT _ q_ F ; q_ C and q ¼ ½T; q ; q T . The specific M; K; U, and F where q_ ¼ T; F C are expressed as follows:
2
3 qc; 0; 0 6 7 M ¼ 4 0; 0; 0 5; 0; 0; s
2
sTt þ Rffqgg n UðqÞ ¼ ; sqt n
$; $
0;
6 K ¼4
F T k$;
I;
ð1 F T Þk$; 0; 2 3 S þ s @S @t 6 7 F¼4 0 5
3
7 0 5; I
ð9Þ
0
where qc; k, and F T are comprised by the material properties, thermal conductivities, and heat conduction model number of each subdomains, respectively; note that the heat source term directly degenerates to S in the Fourier’ s heat conduction since s ¼ 0. B is the derivative of the constraint with respective to the primary variables
B¼
@U @q
ð10Þ
At this point, the first-order time dependent differential algebraic equation has been formulated to describe the heat conduction problem with perfect/imperfect interface. Different heat conduction models can be easily achieved by adjusting the heat conduction model number, F T , and the resulting systems are all in the framework of first-order time systems, which is different from the traditional single-field formulation. It is worth noting that the constraint is only under the expression of the primary variables such that it is classified as Index 2 constraints, which is the highest index existing in the first-order DAE system. 2.3. GS4-1 DAE time integration In this section, the resulting first-order time dependent DAE Index 2 system, Eq. (8), will be solved in the unified time integration framework, termed GS4-1, with the consideration of Index 2 constraint. The specific algorithm will be presented in detail. Let
892
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
qðtÞ : I :¼ ½t 0 ; t f R ! Rndim be the primary unknown vector where t0 and t f are the initial and final time of simulation, respectively. Consider the temporal discretization of the governing equations for a time interval I :¼ ½t 0 ; t f split into sub-intervals, i.e., Sf 1 ½t n ; t nþ1 . The time step size is defined as I :¼ ½t0 ; tf ¼ n¼0 Dtnþ1 :¼ tnþ1 tn > 0, and it is assumed to be constant for simplicity, i.e., Dt ¼ Dtnþ1 . Define the algorithmic governing equation at tnþW 1 ¼ ð1 W 1 Þt n þ W 1 tnþ1 as
~_ þ K f eq e T ~k ¼ Fe ~þB Mq
U qnþ1 ¼ 0
ð11Þ
~_ ¼ q_ n þ W 1 K6 Dq_ q ~ ¼ qn þ W 1 K4 q_ n Dt þ W 2 K5 g1 Dq_ Dt q
ð12Þ
where W i 2 R; Ki 2 Rði ¼ 1; 2; 3; 4; 5; 6Þ and g1 are all algorithmic parameters, which are exploited in the time weighted residual methodology to design the algorithms. More specific derivation can be found in [17]. The updates are given as
qnþ1 ¼ qn þ k4 q_ n Dt þ k5 Dq_ Dt q_ nþ1 ¼ q_ n þ Dq_
ð13Þ
knþ1 ¼ kn þ Dk where ki ði ¼ 1; 2; 3; 4; 5Þ are algorithmic parameters associated with Ki [17]. 2.3.1. Newton-Raphson iteration Employ the Newton-Raphson method to iteratively solve the nonlinear DAE system. At the beginning of each time step, predict the state vectors:
~_ ¼ q_ n þ W 1 K6 Dq_ q ~ ¼ qn þ W 1 K4 q_ n Dt þ W 2 K5 g1 Dq_ Dt q ~k ¼ ð1 W 1 Þkn þ W 1 knþ1
ð14Þ
Define algorithmic Residual vector as
Rg R¼ Rc
"
¼
~_ þ K f eq e T ~k Fe ~þB Mq
U qnþ1
# ¼
0 0
ð15Þ
therefore, the correction can be obtained as kþ1 J Ddkþ1 nþ1 ¼ Rnþ1
q_ nþ1 q_ n Dq_ n ¼ Dk n knþ1 kn
ð16Þ
ð17Þ
and
2 3 eT e þK e J ; W1 B M þ W 2 K5 Dt K W 1 K6 f @R 5 J¼ ¼4
@ Dd L5 DtB qnþ1 ; 0
ð18Þ
e J ¼ @eB T ~k. where K ~ @q Therefore, the correctors are as follows
~_ kþ1 ¼ q ~_ k þ W K Dq_ q 1 6 nþ1 nþ1 nþ1 ~k ~ kþ1 _ q nþ1 ¼ qnþ1 þ W 2 K5 Dqnþ1 Dt ~kkþ1 ¼ ~kk þ W 2 Dknþ1 nþ1
nþ1
Once the solution converges, update the unknowns at the end of the time step is as follow.
~_ kþ1 q_ =W K Dq_ ¼ q n 1 6 nþ1 kþ1
~ Dk ¼ knþ1 kn =W 1 qnþ1 ¼ qn þ k4 q_ n Dt þ k5 Dq_ Dt q_ nþ1 ¼ q_ n þ Dq_
ð21Þ
knþ1 ¼ kn þ Dk 3 þ q11 þ q21 q11 q21 1 ; W 2 K5 ¼ ð1 þ q11 Þð1 þ q21 Þ 2ð1 þ q11 Þð1 þ q21 Þ 1 1 ; k4 ¼ 1; k5 ¼ W 1 ¼ W 1 K4 ¼ ð1 þ q11 Þ 1 þ q21
W 1 K6 ¼
ð22Þ Remark 2.1. 1. The control parameters (q11 and q21 ) are the spurious and principal roots at the high frequency limit which are required to fulfill the following relation
0 6 q11 6 q21 6 1
ð23Þ
2. All the schemes obtained in the framework of GS4-1 are second-order accurate in time, unconditionally stable, zeroorder overshoot, and controllable numerical dissipative [17]. 3. The relationships between the existing algorithms and the GS41 framework are (a) q11 ¼ q21 ¼ 1 represents Crank-Nicolson method; (b) q11 ¼ q21 ¼ 0 represents Gear’s method; (c) q11 ¼ q21 represents the existing algorithms (generalized-a method for first-order transient system) without selective control feature [18]; (d) q11 – q21 represents the new algorithms with selective control feature, which provides numerous schemes with selective dissipation [17]. 3. Numerical examples
where the superscript k in this section represents the numerical iteration step in the Newton-Raphson Iteration and
Dd ¼
ð20Þ
where the specific algorithmic parameters are as follows:
~k ¼ ð1 W 1 Þkn þ W 1 knþ1
until the solution converges
kþ1 kþ1 Ddnþ1 ¼ dnþ1 dknþ1 < tolerance
ð19Þ
In this section, several numerical examples are exploited to verify the effectiveness and illustrate the performance of the proposed method. The numerical examples are all undergone in the following two slabs (1–2) in contact. Numerical Example I comprises the thermal contact of Fouriertype model (F 1T ¼ F 2T ¼ 1) in the slabs and both the perfect interface and the imperfect interface between the two slabs are investigated. The numerical results are compared with the analytical solutions. Numerical Example II comprises the thermal contact between the Fourier-type model (F 1T ¼ 1) and the Cattaneo-model (F 2T ¼ 0) in the two slabs and both the perfect and the imperfect interface are investigated. Numerical Example III considers the thermal contact existing in the two slabs governing by two arbitrary thermal flux models, specifically, F 1T ¼ 0:1; F 2T ¼ 0:5 and k1 ¼ 1; k2 ¼ 2 (see Fig. 1). It’s worth noting that the numerical example II considers the thermal contact between two completely different physics, including (i) the first-order time system associating with the Fourier-type thermal flux model; (ii) the second-order time system associating
893
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
0.2 Slab-2
Slab-1 Interface
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
Fig. 1. Two slabs in contact.
with the Cattaneo-type thermal flux from the perspective of the classical single-field formulation. To achieve the numerical simulations of numerical example II, the time integration is a big challenge. In the traditional single-field formulation, the numerical Time:
1.0 Interface: R=0
0.8 Temperature
0.2 0.4 0.6 0.8 1.0
0.6 0.4 0.2 0.0
0.6
Constraint 1E-16
0.0
example II has to be described by two individual governing equations with respect to each slab, consequently, different time integration schemes have to be considered to deal with these two governing equations. Additionally, these two governing equations not only need to be solved simultaneously from the standpoint of capturing the real physics, but also have to fulfill the interface conditions. Besides, to the authors’ knowledge, numerical examples II and III have not been completely investigated by the existing studies regarding the thermal contact. In this work, the governing equation of heat conduction is illustrated in a generalized two-field formulation, which is a first-order time ODE system. Coupling with the interface conditions, which are expressed in terms of algebraic equations, the ODE becomes a typical first-order time DAE system, in such a way, the proposed GS4-1 DAE time integration can fit and overcome those aforementioned numerical challenges.
ΦT
0.4
Φqf
0.2
Φqc
0.0 -0.2 -0.4 -0.6
0.0
0.2
0.4
0.6
x
0.8
1.0
0.0
0.2
0.4
(a)
1.0
3
Time: 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint 1E-16
Temperature
0.8
(b)
1.0 0.8
0.6 Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
0.6
x
0.8
-3
1.0
0.0
0.2
0.4
(c) Time: 0.2 0.4 0.6 0.8 1.0
1.0
Interface: R=1
0.6 0.4 0.2 0.0
4
Constraint 1E-16
Temperature
0.8
(d)
1.0 0.8
0.6 Time
ΦT
3
Φqf
2
Φqc
1 0 -1 -2
0.0
0.2
0.4
x
(e)
0.6
0.8
1.0
0.0
0.2
0.4
0.6 Time
0.8
1.0
(f)
Fig. 2. The comparison of the thermal contact in the F- and F-slabs. (a, b) are the temperature distribution at different time instants and the constraints in the case of perfect interface (R ¼ 0), respectively; (c, d) are the temperature distribution at different time instants and the constraints in the case of imperfect interface (R ¼ 0:5), respectively; (e, f) are the temperature distribution at different time instants and the constraints in the case of imperfect interface (R ¼ 1), respectively.
894
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
0.6 0.4 0.2
Constraint1E-16
0.2 0.4 0.6 0.8 1.0
Interface: R=0
0.8
Temperature
1.5
Time:
1.0
1.0
ΦT
0.5
Φqc
Φqf
0.0 -0.5 -1.0
0.0 0.0
0.2
0.4
x
0.6
0.8
-1.5
1.0
0.0
0.2
0.4
(a) F 1T
¼ 1 and
F 2T
1.0
¼ 0. (a) the temperature distribution and (b) constraints.
3 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint1E-16
Time:
1.0
Temperature
0.8
(b)
Fig. 3. The perfect thermal contact in the case with
0.8
0.6
Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
x
0.6
0.8
-3
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Time
(a)
(b)
Fig. 4. The imperfect thermal contact (R = 0.5) in the case with F 1T ¼ 1 and F 2T ¼ 0. (a) the temperature distribution and (b) constraints.
The initial and the boundary conditions are given as follows:
Tðx; t ¼ 0Þ ¼ 0; Tð0; tÞ ¼ 1;
Tð1; tÞ ¼ 0
ð24Þ
The GS4-1 DAE integration with parameters (0,0) is exploited in each case. The maximal numerical dissipation is introduced in the time integration. In each simulation, the time interval and the total simulation period are taken as 1E-2 and 1, respectively. Additionally, the consistent Moving Particle Method [19] is directly exploited to achieve the approximations of the different operators (r and r) in Eq. (8), which is given in a strong-form formulation. Due to this procedure, the numerical simulation becomes a relatively simple and direct compared with the weak-form based methods. 3.1. Numerical example I : F 1T ¼ F 2T ¼ 1 In the numerical example I, the thermal contact between two Fourier-type flux governed slabs are investigated. The material and the thermal parameters are taken as q ¼ 1; c ¼ 1, and k1 ¼ k2 ¼ 1. The following cases are treated. 1. R ¼ 0 Thermal contact with perfect interface in which the temperature and the thermal flux are required to be continuous. 2. R ¼ 0:5 Thermal contact with imperfect interface in which the temperature are decided by Eq. (1) while the thermal flux is required to be continuous.
3. R ¼ 1 Thermal contact with imperfect interface. The interface condition is the same as the case of R ¼ 0:5 while a different thermal resistance is considered. As illustrated in Fig. 2, the numerical results (solid lines) of different thermal resistances agree with the analytical solutions (symbols). The constraints with respect to T; qf , and qc are shown in Fig. 2 as well. The continuities of the temperature and flux at the perfect interface case are constrained in the machine error (1E-16). Similarly, the fulfillment of the constraints in the imperfect interface cases also show good performance. Due to the difference of the thermal resistance, the temperature distribution illustrates different physics. As evident from the numerical results, the larger thermal resistance leads to a large temperature jump at the interface, which agrees with the practical physics. 3.2. Numerical example II : F 1T ¼ 1; F 2T ¼ 0 In this part, we consider a very special thermal contact case between the Fourier-type (F T ¼ 1) heat conduction and the Cattaneo-type (F T ¼ 0) heat conduction. This case usually happens in the engineering problems when the thermal contact is considered in multi-layered films and the thermal interaction between film and its subtract. In the classical single-field formulation, the two slabs have to be modeled separately and different kinds of time integrations also need to be considered. The resulting two governing equations are fundamentally different due to the difference between the parabolic heat conduction and the hyperbolic
895
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
Time: Interface: R=0
0.8
Temperature
0.2 0.4 0.6 0.8 1.0
0.6 0.4 0.2
0.6
Constraint1E-16
1.0
ΦT
0.4
Φqf
0.2
Φqc
0.0 -0.2 -0.4
0.0
-0.6 0.0
0.2
0.4
0.6
x
0.8
1.0
0.0
0.2
0.4
(a)
1.0
3
Time: 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint1E-16
Temperature
0.8
(b)
1.0 0.8
0.6
Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
x
0.6
0.8
1.0
-3
0.0
0.2
0.4
0.6
0.8
1.0
Time
(c)
(d)
Fig. 5. The imperfect thermal contact (R = 0.5) in the case with F 1T ¼ 0:1 and F 2T ¼ 0:5. (a) the temperature distribution and (b) constraints.
heat conduction. Additionally, the thermal flux of the hyperbolic heat condition is given implicitly with respect to the temperature in the single-field formulation such that it is difficult to obtain the constraint equation in the framework of the DAE system. Unlike the traditional single-field form, both the thermal flux and temperature in the proposed model can be easily obtained and the associated constraint equations are also formulated in a concise way. Fig. 3 illustrates the perfect thermal contact between the Fourier-type heat conduction and Cattaneo-type heat conduction. Different heat propagation patterns can be observed in the two slabs, specifically, the heat signal propagates with an infinite speed such that the temperature can be observed at any place of Slab-1; while the heat signal propagates with an finite speed in the Slab-2 such that the telegraphic propagation can be observed. Fig. 4 presents the case with imperfect thermal interface (R ¼ 0:5). Similarly, two different types of heat conductions can be observed in the temperature distribution at each time instants. However, compared with Fig. 3, the temperature of the Fourier’ s slab increases faster while the temperature of the Cattaneo’s slab increases slower than that of the perfect thermal contact case. This is because the heat transfer from the Fourier’s region to the Cattaneo’s region is resisted by the imperfect interface in such a way more energy in the Fourier’s region is reserved to increase the temperature of the Slab 1 instead of transferring to the Slab 2. Besides, the constraints are excellent preserved in the machine error (1E-16). 3.3. Numerical example III : F 1T ¼ 0:1; F 2T ¼ 0:5 In this case, we consider the thermal contact of a general case emanating from the C-F heat conduction model. The heat conduc-
tivities are also considered to be different of the two slabs as: k1 ¼ 1 and k2 ¼ 2. Different heat conduction processes can be observed in Fig. 5. In general, the numerical results reveal the underlying physics of different F T . It’s worth noting that the parameter F T represents the portion of the Fourier-type in such a way the feature of telegraphic propagation in this case is stronger in the Slab-1 with smaller F T . In addition, Fig. 5(b) in the perfect thermal contact case shows that the constraints are automatically achieved at the beginning of the simulation. This may be due to that the heat signal is very weak at the interface such that the constraint does not play a role in the algorithm. Also, most of the heat is stored in the Slab-1 due to the thermal resistance, which leads the temperature of Slab-1 increases drastically while the heat of Slab-2 illustrates a gentle transfer. 4. Concluding remarks In this work, we have presented the formalism of a generalized heat conduction model with the consideration of perfect/imperfect thermal interface, which unifies both the Cattaneo’s and Jeffreys’ heat conduction (second-order time dependent system in terms of the single-field form) and the Fourier ’s type heat conduction (first order time dependent system in terms of the single-field form) into a simple two-field first-order time dependent differential algebraic framework in which the interface conditions are treated as algebraic equations, classified as Index 2 Constraints. The GS4 DAE Index 2 time integration is applied to carry out the numerical simulations. The numerical results from the GS4 DAE Index 2 time integration show the particular features of the thermal contacts between either Fourier- and/or Cattaneo- and/or Jeffreys- types with perfect/imperfect thermal interfaces. The proposed method has be proved to have the ability of describing
896
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
the real heat conduction existing in the engineering problems such as multi-layer composites, film-substrate structures, and bio-tissue system, etc. A direct application of the proposed computational framework can be taken to is that the heat transport in the skins, which are usually comprised by epidermis, dermis, and hypodemis [20,21]. Each layer has its own material properties and relaxation times, which is a practical problem of the numerical example III. To model the heat process in such complex problem, the proposed computational framework can be directly implemented. Moreover, the structure of a general optical device is usually comprised by thin optical films and a thick subtract. The heat conduction of the film-subtract system has been described by the heat conduction model with relaxation time in the previous studies [9]. As a matter of fact, the thickness of the thin film may be at a close order of magnitude of its heat carrier’s mean free path, whereas the size of the subtract is usually much larger than the mean free path of its heat carriers, especially for thick metal subtracts, such as Fe and Cu. Using the same-type of heat conduction model seems to lack the reliability. The coupling different types of heat conduction models, such as the thermal contact between Fourier’ s and non-Fourier’ s heat conduction, appears to be a computational challenge in the previous studies since the heat flux continuity is difficult to be achieved. However, it has been proved to be very easy for the proposed computational framework and the numerical examples II and III have revealed this particular feature of the proposed method. The further works will mainly focus on the following two aspects: (i) the thermal contact in highly heterogeneous systems and the cases which are small in both the space and time; (ii) the application of general interface conditions on various imperfect thermal contacts. Acknowledgments We wish to thank the reviewers for helpful suggestions. We also would like to thank Prof. Kumar K. Tamma, Dr. Masao Shimada, and Allon Zhu from the University of Minnesota for their advice and support. Conflict of interest Authors declare that there is no conflict of interest. References [1] D. Hasselman, L.F. Johnson, Effective thermal conductivity of composites with interfacial thermal barrier resistance, J. Compos. Mater. 21 (6) (1987) 508–515.
[2] S. Torquato, M. Rintoul, Effect of the interface on the properties of composite media, Phys. Rev. Lett. 75 (22) (1995) 4067. [3] J. Yvonnet, Q.-C. He, Q.-Z. Zhu, J.-F. Shao, A general and efficient computational procedure for modelling the kapitza thermal resistance based on xfem, Comput. Mater. Sci. 50 (4) (2011) 1220–1224. [4] A. Javili, S. Kaessmair, P. Steinmann, General imperfect interfaces, Comput. Methods Appl. Mech. Eng. 275 (2014) 76–97. [5] A. Jain, K. Tamma, Parabolic heat conduction specialized applications involving imperfect contact surfaces: local discontinuous galerkin finite element method-part 2, J. Therm. Stresses 33 (4) (2010) 344–355. [6] H.-T. Chen, K.-C. Liu, Study of hyperbolic heat conduction problem in the film and substrate composite with the interface resistance, Jpn. J. Appl. Phys. 41 (10R) (2002) 6267. [7] D.Y. Tzou, Reflection and refraction of thermal waves from a surface or an interface between dissimilar materials, Int. J. Heat Mass Transfer 36 (2) (1993) 401–410. [8] K.-C. Liu, Analysis of dual-phase-lag thermal behaviour in layered films with temperature-dependent interface thermal resistance, J. Phys. D: Appl. Phys. 38 (19) (2005) 3722. [9] K.-C. Liu, Analysis of thermal behavior in multi-layer metal thin films based on hyperbolic two-step model, Int. J. Heat Mass Transfer 50 (7) (2007) 1397– 1407. [10] 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. [11] K.-C. Liu, H.-T. Chen, Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment, Int. J. Heat Mass Transfer 52 (5) (2009) 1185–1192. [12] G. Yeoh, X. Gu, V. Timchenko, S. Valenzuela, B. Cornell, High order accurate dual-phase-lag numerical model for microscopic heating in multiple domains, Int. Commun. Heat Mass Transfer 78 (2016) 21–28. [13] W. Little, The transport of heat between dissimilar solids at low temperatures, Can. J. Phys. 37 (3) (1959) 334–349. [14] K.K. Tamma, X. Zhou, Macroscale and microscale thermal transport and thermo-mechanical interactions: some noteworthy perspectives, J. Therm. Stresses 21 (3–4) (1998) 405–449. [15] X. Zhou, K. Tamma, A new unified theory underlying time dependent linear first-order systems: a prelude to algorithms by design, Int. J. Numer. Meth. Eng. 60 (10) (2004) 1699–1740. [16] K. Tamma, D. Sha, X. Zhou, Time discretized operators. Part 1: Towards the theoretical design of a new generation of a generalized family of unconditionally stable implicit and explicit representations of arbitrary order for computational dynamics, Comput. Meth. Appl. Mech. Eng. 192 (3) (2003) 257–290. [17] S. Masuri, M. Sellier, X. Zhou, K. Tamma, Design of order-preserving algorithms for transient first-order systems with controllable numerical dissipation, Int. J. Numer. Meth. Eng. 88 (13) (2011) 1411–1448. [18] K.E. Jansen, C.H. Whiting, G.M. Hulbert, A generalized-a method for integrating the filtered navier–stokes equations with a stabilized finite element method, Comput. Meth. Appl. Mech. Eng. 190 (3) (2000) 305–319. [19] T. Xue, K.K. Tamma, X. Zhang, A consistent moving particle system simulation method: applications to parabolic/hyperbolic heat conduction type problems, Int. J. Heat Mass Transfer 101 (2016) 365–372. [20] F. Xu, K. Seffen, T. Lu, Non-fourier analysis of skin biothermomechanics, Int. J. Heat Mass Transfer 51 (9) (2008) 2237–2259. [21] J.T. Whitton, J. Everall, The thickness of the epidermis, Br. J. Dermatol. 89 (5) (1973) 467–476.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Generalized heat conduction model involving imperfect thermal contact surface: Application of the GSSSS-1 differential-algebraic equation time integration Tao Xue a, Xiaobing Zhang a,⇑, Kumar K. Tamma b 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, Minneapolis, USA
a r t i c l e
i n f o
Article history: Received 6 May 2017 Received in revised form 19 September 2017 Accepted 21 September 2017 Available online 17 October 2017 Keywords: Thermal contact Imperfect interface Differential algebraic equation Generalized heat conduction model
a b s t r a c t A number of contributions have been made to model thermal contact problems involving perfect/imperfect thermal interfaces. The Fourier-, Cattaneo-, and Jeffreys-types of thermal flux modes have been exploited to govern the heat transfer processes of adjacent regions. However, most of the existing studies consider using only one type of these thermal flux models to describe the thermal physics of the entire system such that the mathematical model may fail to precisely describe the physics within the multidomain system, especially, when the mean free paths of different adjacent domains span in a large range with respect to the characteristic dimension of the material. To circumvent this issue, a generalized heat conduction model, termed C- and F-processes heat conduction model, with the consideration of perfect/ imperfect thermal interface is proposed and formulated under the expression of a first-order time dependent differential-algebraic equation in which the interface conditions are treated as algebraic equations. The unified time integration of GSSSS-1 is extended to solve the resulting differential-algebraic system with Index 2 constraints. The numerical results illustrate that the proposed framework has a better capacity of describing the thermal contacts with/without the thermal resistance of Fourier-Fourier, Fourier-Cattaneo, and other general cases. Ó 2017 Published by Elsevier Ltd.
1. Introduction For the coupling of different scales/materials within a solid, of either multilayered nano-scale films, micro devices or macro systems, the interfaces connecting adjacent domains play an important role in the thermal response of the entire system. In the ideal situation, the adjacent domains contact perfectly such that the temperature and the flux perform no jumps across the interface due to the energy conservation law. The interface of the perfect contacting case is classified as Perfect Interface. However, in the general practical applications, the adjacent domains usually do not contact well and there exist gaps between the contacting domains. In such a way, this Imperfect Interface leads a discontinuity of the temperature while the thermal flux retains a continuity across the interface. Specifically, the relationships with respect to the temperature T and the thermal flux q for the imperfect interfaces [1–3] are given by: ⇑ Corresponding author. E-mail addresses: [email protected] (T. Xue), [email protected] (X. Zhang), [email protected] (K.K. Tamma). https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.081 0017-9310/Ó 2017 Published by Elsevier Ltd.
sTt ¼ Rffqgg n sqt n ¼ 0
ð1Þ
where sTt represents TjC1 TjC2 and ffqgg represents qjC1 þ qjC2 . The coefficient R P 0 represents the thermal (Kapitza)
1 2
resistance, when R ¼ 0 the thermal contact is associated with the perfect interface while R – 0 the thermal contact is associated with the imperfect interface. C1 and C2 are the surfaces in contact and n is the normal direction from C2 to C1 . It is worth noting that Eq. (1) represents the lowly-conducting interface, which is classically modeled using Kapitza’s assumption of thermal resistance. A more general formulation regarding the interface conditions has been proposed in [4], both the lowly-conducting and the highlyconducting interfaces are covered in this general imperfect interface theory. However, the general imperfect interface assumes that the interface is a zero-thickness material and has its own heat transfer process such that it may introduce more numerical problems. For simplicity, the lowly-conducting interface is of interest to this article. Numerous analytical, experimental, and numerical models have been done regarding the thermal contact and the treatments of the
890
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
interface condition. In some thermal contact problems, the adjacent regions are governed by the Fourier’ s heat conduction model. The weak-form based numerical approaches such as FEM, XFEM are applied to capture the strong discontinuity caused by the thermal resistance [5,1,3]. In the engineering applications of thermal wave theory involving thermal interfaces, the thermal contact is usually considered in the framework of the hyperbolic heat conduction model. The thermal resistance was treated as a jump boundary conditions in [6,7]. Additionally, remarkable studies regarding the theoretical analysis with respect to the thermal contact with/without resistance in the framework of hyperbolic model have been done in [8,9]. Note that the numerical approaches applied in the hyperbolic-type system should be able to capture the discontinuities in each domain and also works well under the condition of thermal resistance. The thermal contact has also been extended to the field of bio-engineering, such as the thermal processes in the multi-layer bio-tissues. Both the Fourier-type heat conduction and dual-phase-lag heat conduction model1 have been done with respect to this particular aspect [11,12]. Instead of the linear thermal resistance between the temperature and the flux at the interface, the heat flux across the interface of two layers is assumed to proportionally change with the difference between the fourth powers of the temperatures on the two sides of the interface, which is similar to the radiation heat transfer model. The specific expression of R [8] is as follows:
TjC TjC2 R¼ 1 T4j T4j h C1 C2
ð2Þ
is a thermal contact parameter related to the surface where h roughness [13]. The overview described earlier indicates that the existing works consider implementing only one certain thermal flux model to govern the thermal transfer over the entire multi-domain system, which seems to be not precise and cannot be considered as a complete theory when the mean free paths of different adjacent domains span in a large range with respect to the characteristic dimension of the material. For example, [6] utilized the hyperbolic heat conduction model in all the adjacent regions to formulate the thermal contact of the practical electronic device, which are comprised by conducting film and substrate composite. However, the proper thermal flux models for the film and the subtract are necessary to be determined by the relationship between the mean free path (k) and the characteristic dimension of the material (L). Specifically, 1. When k L, the pure ballistic transport dominates the transport process; 2. When k L, the transport process is governed by a diffusiveballistic nature; 3. When k L, the heat transport is treat to be macroscopic and is governed by a purely diffusive nature. In other words, using a single heat conduction model cannot fulfill the requirement of precisely describing the thermal process over the entire system. In order to achieve it, different types of heat conduction models are coupled by the interface, which turns out to be a numerical challenge for the previous studies. The specific reasons include that: (i) different orders of time-dependent systems are needed to describe the thermal contact between the Fourier-type flux and the Cattaneo-/Jeffreys-type flux; 1 The underlying physics of the dual-phase-lag has been proved to be the same as the Jeffreys’ conduction model. More details can be found in [10].
(ii) the flux models of the Cattaneo-type and the Jeffreys-type involve the partial time derivative and are not explicitly given, which implies that the interface condition of thermal flux is not easy to be fulfilled; (iii) the existing numerical treatments with respect to both the perfect and imperfect interfaces also seem to be a numerical challenge in different approaches, specifically, the interface is usually handled as a boundary condition that are embedded in the time stepping of the governing equation. Consequently, in most of the previous studies the governing equations of each domain have to be evaluated separately and the time stepping strategies to the associated numerical simulations are designed in the framework of staggering methodology, which may have the issue of violating the interface conditions. Therefore, increasing application and necessity regarding the treatment of interface, the modeling of the heat conduction within adjacent regions, and the thermal performance at different scales have motivated the objective of this work to develop a general framework of describing the complex thermal contact physics and a decent approach to deal with the perfect/imperfect interface condition. To circumvent these issues, a two-field (both the temperature and the thermal flux are primary variables) generalized heat conduction model, the C- and F-processes heat conduction model, is utilized to depict the heat conduction of the adjacent regions. The C- and F-processes heat conduction model is originally developed in [14] and is based upon the assumption that the Cattaneo-type thermal flux associated with a finite sound speed and the Fourier-type thermal flux associated with an infinite sound speed co-exist in a heat process simultaneously. Therefore, the C- and F-heat conduction model is acknowledged to cover various heat conduction models. Though the hyperbolicand Jeffreys-types of heat conduction are acknowledged as second-order time systems in terms of single-field formulation (only the temperature is the primary variable), the governing equation of the proposed C-F heat conduction is presented in a first-order two-field form with considering both the temperature and the thermal flux as the primary physical fields of interest. In the framework of the two-field C-and F-model, the thermal flux model in each adjacent region can be switched easily and all the governing equations can be expressed in a unified formulation. Different from the past efforts, we treat the interface condition, Eq. (1), as an algebraic equation such that coupling the interface algebraic equation with the generalized two-field C- and F-heat conduction equation formulates a first-order time dependent Differential Algebraic Equation. The principal contributions emanating from such a generalized framework are summarized as follows: 1. A simple computational framework to model the thermal contact between either the Fourier-type and/or Cattaneo-type and/or Jeffreys-type, without having to resort to the individual governing equation. 2. A perfect preservation of the perfect/imperfect thermal interface condition in a machine error (1E-16). Additionally, we also present the formalism of a unified secondorder accurate time integration, GS4 (generalized single step single solve)[15–17], to solve this first order DAEs system in terms of GS4-1 DAE Index-2. Note that the Index 2 is the highest constraint index existing in the first-order time dependent DAE system and naturally arises from the interface condition in this particular thermal contact problem.
891
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
2. F T 2 ð0; 1Þ: Jeffreys-type Model
2. Theory In this section, we first introduce the fundamental theory of the generalized C- and F-processes model and its associated two-field form. Then, the differential algebraic equation is formulated to describe the heat conduction model involving the perfect/imperfect thermal interface. Last, the GS4-1 DAE Index 2 time integration is derived in detail.
The standard C-and F-heat conduction model was developed from the fundamental Boltzmann Transport Equation under the assumption that the Fourier-type and the Cattaneo-type processes coexist simultaneously in the heat conduction process [14,10]. Compared with the Jeffreys-type heat conduction model (DualPhase-Lag heat conduction), which can only reduce to a Fourierlike heat conduction model with relaxation, the C- and F-model not only explains the underlying theory of the slow processes (Cattaneo-type), fast processes (Fourier-type), and those in between; but also naturally leads to a generalization of the heat conduction in solids of the Jeffreys-, Cattaneo-, and Fourier-type models. The standard C-and F-model is given by:
ð3Þ
where s represents the relaxation time; q; qF , and qC are the total heat flux, the fast heat flux (Fourier-type), and the slow heat flux (Cattaneo-type), respectively. The corresponding heat conductiviF ties of q; qF , and qC are k; kF , and kC , respectively. F T ¼ kCKþk is defined F as a macro-scale heat conduction model number to physically illustrate different types of heat conduction processes, specifically, 1. F T ¼ 0 The C- and F-model degenerates to the Cattaneo-type (hyperbolic) heat conduction processes, which has a finite heat propagation speed. 2. F T ¼ 1 The C- and F-model degenerates to the Fourier-type (parabolic) heat conduction processes, which has an infinite heat propagation speed. 3. F T 2 ð0; 1Þ The C- and F-model degenerates to the Jeffreys-type (parabolic) heat conduction processes, the discontinuities existing in the Jeffreys-type model is smoothed by the diffusion (Fouriertype) effects. It’s worth noting that the acknowledged Jeffreystype heat conduction is also considered as a generalized model, however, it’s not able to obtain the Fourier-type heat conduction. The conservation of energy yields the following
1 @T ¼rqþS a @t
ð4Þ
where a ¼ qkc ; q and c represent the material density and specific heat, respectively; S is a heat source within the material. Substituting the C-and F-heat flux model to Eq. (4), the corresponding generalized single-field form with respect to the one temperature theory are illustrated as follows: 1. F T ¼ 0: Cattaneo-type Model
1 @ T 1 @T 1 @S 2 þ r T þ S þ s ¼ @t s c2T @t2 a @t 2
ð5Þ
ð6Þ
3. F T ¼ 1: Fourier-type Model
1 @T
a @t
¼ r2 T þ S
where cT ¼
2.1. Generalized heat conduction C- and F-model
qF ¼ kF T rT @q qC þ s C ¼ ð1 F T ÞkrT @t q ¼ qF þ qC
1 @ 2 T 1 @T @ r2 T 1 @S 2 þ r T þ s F S þ s ¼ þ T @t @t s c2T @t 2 a @t
paffiffi s
ð7Þ
in the Cattaneo- and Jeffreys-type models.
2.2. Differential algebraic equations of thermal contact based on the Cand F-heat conduction model In this work, instead of the traditional single-field form, we propose to directly utilize the two-field form of the C- and F-model in the aspect of the thermal contact problem. Both the temperature and the thermal flux are treated as the primary variables. The interface conditions are naturally comprised by these two primary variables and considered to be the constraints with respect to the governing equation. To formulate the DAE with Index 2 constraint, the interface conditions are treated as constraint and a Lagrangian Multiplier (k) is introduced to couple the constraints with the generalized C- and F-heat conduction model. The specific formulation is as follows:
M q_ þ Kq þ BT k ¼ F UðqÞ ¼ 0
ð8Þ
h iT _ q_ F ; q_ C and q ¼ ½T; q ; q T . The specific M; K; U, and F where q_ ¼ T; F C are expressed as follows:
2
3 qc; 0; 0 6 7 M ¼ 4 0; 0; 0 5; 0; 0; s
2
sTt þ Rffqgg n UðqÞ ¼ ; sqt n
$; $
0;
6 K ¼4
F T k$;
I;
ð1 F T Þk$; 0; 2 3 S þ s @S @t 6 7 F¼4 0 5
3
7 0 5; I
ð9Þ
0
where qc; k, and F T are comprised by the material properties, thermal conductivities, and heat conduction model number of each subdomains, respectively; note that the heat source term directly degenerates to S in the Fourier’ s heat conduction since s ¼ 0. B is the derivative of the constraint with respective to the primary variables
B¼
@U @q
ð10Þ
At this point, the first-order time dependent differential algebraic equation has been formulated to describe the heat conduction problem with perfect/imperfect interface. Different heat conduction models can be easily achieved by adjusting the heat conduction model number, F T , and the resulting systems are all in the framework of first-order time systems, which is different from the traditional single-field formulation. It is worth noting that the constraint is only under the expression of the primary variables such that it is classified as Index 2 constraints, which is the highest index existing in the first-order DAE system. 2.3. GS4-1 DAE time integration In this section, the resulting first-order time dependent DAE Index 2 system, Eq. (8), will be solved in the unified time integration framework, termed GS4-1, with the consideration of Index 2 constraint. The specific algorithm will be presented in detail. Let
892
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
qðtÞ : I :¼ ½t 0 ; t f R ! Rndim be the primary unknown vector where t0 and t f are the initial and final time of simulation, respectively. Consider the temporal discretization of the governing equations for a time interval I :¼ ½t 0 ; t f split into sub-intervals, i.e., Sf 1 ½t n ; t nþ1 . The time step size is defined as I :¼ ½t0 ; tf ¼ n¼0 Dtnþ1 :¼ tnþ1 tn > 0, and it is assumed to be constant for simplicity, i.e., Dt ¼ Dtnþ1 . Define the algorithmic governing equation at tnþW 1 ¼ ð1 W 1 Þt n þ W 1 tnþ1 as
~_ þ K f eq e T ~k ¼ Fe ~þB Mq
U qnþ1 ¼ 0
ð11Þ
~_ ¼ q_ n þ W 1 K6 Dq_ q ~ ¼ qn þ W 1 K4 q_ n Dt þ W 2 K5 g1 Dq_ Dt q
ð12Þ
where W i 2 R; Ki 2 Rði ¼ 1; 2; 3; 4; 5; 6Þ and g1 are all algorithmic parameters, which are exploited in the time weighted residual methodology to design the algorithms. More specific derivation can be found in [17]. The updates are given as
qnþ1 ¼ qn þ k4 q_ n Dt þ k5 Dq_ Dt q_ nþ1 ¼ q_ n þ Dq_
ð13Þ
knþ1 ¼ kn þ Dk where ki ði ¼ 1; 2; 3; 4; 5Þ are algorithmic parameters associated with Ki [17]. 2.3.1. Newton-Raphson iteration Employ the Newton-Raphson method to iteratively solve the nonlinear DAE system. At the beginning of each time step, predict the state vectors:
~_ ¼ q_ n þ W 1 K6 Dq_ q ~ ¼ qn þ W 1 K4 q_ n Dt þ W 2 K5 g1 Dq_ Dt q ~k ¼ ð1 W 1 Þkn þ W 1 knþ1
ð14Þ
Define algorithmic Residual vector as
Rg R¼ Rc
"
¼
~_ þ K f eq e T ~k Fe ~þB Mq
U qnþ1
# ¼
0 0
ð15Þ
therefore, the correction can be obtained as kþ1 J Ddkþ1 nþ1 ¼ Rnþ1
q_ nþ1 q_ n Dq_ n ¼ Dk n knþ1 kn
ð16Þ
ð17Þ
and
2 3 eT e þK e J ; W1 B M þ W 2 K5 Dt K W 1 K6 f @R 5 J¼ ¼4
@ Dd L5 DtB qnþ1 ; 0
ð18Þ
e J ¼ @eB T ~k. where K ~ @q Therefore, the correctors are as follows
~_ kþ1 ¼ q ~_ k þ W K Dq_ q 1 6 nþ1 nþ1 nþ1 ~k ~ kþ1 _ q nþ1 ¼ qnþ1 þ W 2 K5 Dqnþ1 Dt ~kkþ1 ¼ ~kk þ W 2 Dknþ1 nþ1
nþ1
Once the solution converges, update the unknowns at the end of the time step is as follow.
~_ kþ1 q_ =W K Dq_ ¼ q n 1 6 nþ1 kþ1
~ Dk ¼ knþ1 kn =W 1 qnþ1 ¼ qn þ k4 q_ n Dt þ k5 Dq_ Dt q_ nþ1 ¼ q_ n þ Dq_
ð21Þ
knþ1 ¼ kn þ Dk 3 þ q11 þ q21 q11 q21 1 ; W 2 K5 ¼ ð1 þ q11 Þð1 þ q21 Þ 2ð1 þ q11 Þð1 þ q21 Þ 1 1 ; k4 ¼ 1; k5 ¼ W 1 ¼ W 1 K4 ¼ ð1 þ q11 Þ 1 þ q21
W 1 K6 ¼
ð22Þ Remark 2.1. 1. The control parameters (q11 and q21 ) are the spurious and principal roots at the high frequency limit which are required to fulfill the following relation
0 6 q11 6 q21 6 1
ð23Þ
2. All the schemes obtained in the framework of GS4-1 are second-order accurate in time, unconditionally stable, zeroorder overshoot, and controllable numerical dissipative [17]. 3. The relationships between the existing algorithms and the GS41 framework are (a) q11 ¼ q21 ¼ 1 represents Crank-Nicolson method; (b) q11 ¼ q21 ¼ 0 represents Gear’s method; (c) q11 ¼ q21 represents the existing algorithms (generalized-a method for first-order transient system) without selective control feature [18]; (d) q11 – q21 represents the new algorithms with selective control feature, which provides numerous schemes with selective dissipation [17]. 3. Numerical examples
where the superscript k in this section represents the numerical iteration step in the Newton-Raphson Iteration and
Dd ¼
ð20Þ
where the specific algorithmic parameters are as follows:
~k ¼ ð1 W 1 Þkn þ W 1 knþ1
until the solution converges
kþ1 kþ1 Ddnþ1 ¼ dnþ1 dknþ1 < tolerance
ð19Þ
In this section, several numerical examples are exploited to verify the effectiveness and illustrate the performance of the proposed method. The numerical examples are all undergone in the following two slabs (1–2) in contact. Numerical Example I comprises the thermal contact of Fouriertype model (F 1T ¼ F 2T ¼ 1) in the slabs and both the perfect interface and the imperfect interface between the two slabs are investigated. The numerical results are compared with the analytical solutions. Numerical Example II comprises the thermal contact between the Fourier-type model (F 1T ¼ 1) and the Cattaneo-model (F 2T ¼ 0) in the two slabs and both the perfect and the imperfect interface are investigated. Numerical Example III considers the thermal contact existing in the two slabs governing by two arbitrary thermal flux models, specifically, F 1T ¼ 0:1; F 2T ¼ 0:5 and k1 ¼ 1; k2 ¼ 2 (see Fig. 1). It’s worth noting that the numerical example II considers the thermal contact between two completely different physics, including (i) the first-order time system associating with the Fourier-type thermal flux model; (ii) the second-order time system associating
893
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
0.2 Slab-2
Slab-1 Interface
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
Fig. 1. Two slabs in contact.
with the Cattaneo-type thermal flux from the perspective of the classical single-field formulation. To achieve the numerical simulations of numerical example II, the time integration is a big challenge. In the traditional single-field formulation, the numerical Time:
1.0 Interface: R=0
0.8 Temperature
0.2 0.4 0.6 0.8 1.0
0.6 0.4 0.2 0.0
0.6
Constraint 1E-16
0.0
example II has to be described by two individual governing equations with respect to each slab, consequently, different time integration schemes have to be considered to deal with these two governing equations. Additionally, these two governing equations not only need to be solved simultaneously from the standpoint of capturing the real physics, but also have to fulfill the interface conditions. Besides, to the authors’ knowledge, numerical examples II and III have not been completely investigated by the existing studies regarding the thermal contact. In this work, the governing equation of heat conduction is illustrated in a generalized two-field formulation, which is a first-order time ODE system. Coupling with the interface conditions, which are expressed in terms of algebraic equations, the ODE becomes a typical first-order time DAE system, in such a way, the proposed GS4-1 DAE time integration can fit and overcome those aforementioned numerical challenges.
ΦT
0.4
Φqf
0.2
Φqc
0.0 -0.2 -0.4 -0.6
0.0
0.2
0.4
0.6
x
0.8
1.0
0.0
0.2
0.4
(a)
1.0
3
Time: 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint 1E-16
Temperature
0.8
(b)
1.0 0.8
0.6 Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
0.6
x
0.8
-3
1.0
0.0
0.2
0.4
(c) Time: 0.2 0.4 0.6 0.8 1.0
1.0
Interface: R=1
0.6 0.4 0.2 0.0
4
Constraint 1E-16
Temperature
0.8
(d)
1.0 0.8
0.6 Time
ΦT
3
Φqf
2
Φqc
1 0 -1 -2
0.0
0.2
0.4
x
(e)
0.6
0.8
1.0
0.0
0.2
0.4
0.6 Time
0.8
1.0
(f)
Fig. 2. The comparison of the thermal contact in the F- and F-slabs. (a, b) are the temperature distribution at different time instants and the constraints in the case of perfect interface (R ¼ 0), respectively; (c, d) are the temperature distribution at different time instants and the constraints in the case of imperfect interface (R ¼ 0:5), respectively; (e, f) are the temperature distribution at different time instants and the constraints in the case of imperfect interface (R ¼ 1), respectively.
894
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
0.6 0.4 0.2
Constraint1E-16
0.2 0.4 0.6 0.8 1.0
Interface: R=0
0.8
Temperature
1.5
Time:
1.0
1.0
ΦT
0.5
Φqc
Φqf
0.0 -0.5 -1.0
0.0 0.0
0.2
0.4
x
0.6
0.8
-1.5
1.0
0.0
0.2
0.4
(a) F 1T
¼ 1 and
F 2T
1.0
¼ 0. (a) the temperature distribution and (b) constraints.
3 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint1E-16
Time:
1.0
Temperature
0.8
(b)
Fig. 3. The perfect thermal contact in the case with
0.8
0.6
Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
x
0.6
0.8
-3
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Time
(a)
(b)
Fig. 4. The imperfect thermal contact (R = 0.5) in the case with F 1T ¼ 1 and F 2T ¼ 0. (a) the temperature distribution and (b) constraints.
The initial and the boundary conditions are given as follows:
Tðx; t ¼ 0Þ ¼ 0; Tð0; tÞ ¼ 1;
Tð1; tÞ ¼ 0
ð24Þ
The GS4-1 DAE integration with parameters (0,0) is exploited in each case. The maximal numerical dissipation is introduced in the time integration. In each simulation, the time interval and the total simulation period are taken as 1E-2 and 1, respectively. Additionally, the consistent Moving Particle Method [19] is directly exploited to achieve the approximations of the different operators (r and r) in Eq. (8), which is given in a strong-form formulation. Due to this procedure, the numerical simulation becomes a relatively simple and direct compared with the weak-form based methods. 3.1. Numerical example I : F 1T ¼ F 2T ¼ 1 In the numerical example I, the thermal contact between two Fourier-type flux governed slabs are investigated. The material and the thermal parameters are taken as q ¼ 1; c ¼ 1, and k1 ¼ k2 ¼ 1. The following cases are treated. 1. R ¼ 0 Thermal contact with perfect interface in which the temperature and the thermal flux are required to be continuous. 2. R ¼ 0:5 Thermal contact with imperfect interface in which the temperature are decided by Eq. (1) while the thermal flux is required to be continuous.
3. R ¼ 1 Thermal contact with imperfect interface. The interface condition is the same as the case of R ¼ 0:5 while a different thermal resistance is considered. As illustrated in Fig. 2, the numerical results (solid lines) of different thermal resistances agree with the analytical solutions (symbols). The constraints with respect to T; qf , and qc are shown in Fig. 2 as well. The continuities of the temperature and flux at the perfect interface case are constrained in the machine error (1E-16). Similarly, the fulfillment of the constraints in the imperfect interface cases also show good performance. Due to the difference of the thermal resistance, the temperature distribution illustrates different physics. As evident from the numerical results, the larger thermal resistance leads to a large temperature jump at the interface, which agrees with the practical physics. 3.2. Numerical example II : F 1T ¼ 1; F 2T ¼ 0 In this part, we consider a very special thermal contact case between the Fourier-type (F T ¼ 1) heat conduction and the Cattaneo-type (F T ¼ 0) heat conduction. This case usually happens in the engineering problems when the thermal contact is considered in multi-layered films and the thermal interaction between film and its subtract. In the classical single-field formulation, the two slabs have to be modeled separately and different kinds of time integrations also need to be considered. The resulting two governing equations are fundamentally different due to the difference between the parabolic heat conduction and the hyperbolic
895
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
Time: Interface: R=0
0.8
Temperature
0.2 0.4 0.6 0.8 1.0
0.6 0.4 0.2
0.6
Constraint1E-16
1.0
ΦT
0.4
Φqf
0.2
Φqc
0.0 -0.2 -0.4
0.0
-0.6 0.0
0.2
0.4
0.6
x
0.8
1.0
0.0
0.2
0.4
(a)
1.0
3
Time: 0.2 0.4 0.6 0.8 1.0
Interface: R=0.5
0.6 0.4 0.2
Constraint1E-16
Temperature
0.8
(b)
1.0 0.8
0.6
Time
2
ΦT
1
Φqc
Φqf
0 -1 -2
0.0 0.0
0.2
0.4
x
0.6
0.8
1.0
-3
0.0
0.2
0.4
0.6
0.8
1.0
Time
(c)
(d)
Fig. 5. The imperfect thermal contact (R = 0.5) in the case with F 1T ¼ 0:1 and F 2T ¼ 0:5. (a) the temperature distribution and (b) constraints.
heat conduction. Additionally, the thermal flux of the hyperbolic heat condition is given implicitly with respect to the temperature in the single-field formulation such that it is difficult to obtain the constraint equation in the framework of the DAE system. Unlike the traditional single-field form, both the thermal flux and temperature in the proposed model can be easily obtained and the associated constraint equations are also formulated in a concise way. Fig. 3 illustrates the perfect thermal contact between the Fourier-type heat conduction and Cattaneo-type heat conduction. Different heat propagation patterns can be observed in the two slabs, specifically, the heat signal propagates with an infinite speed such that the temperature can be observed at any place of Slab-1; while the heat signal propagates with an finite speed in the Slab-2 such that the telegraphic propagation can be observed. Fig. 4 presents the case with imperfect thermal interface (R ¼ 0:5). Similarly, two different types of heat conductions can be observed in the temperature distribution at each time instants. However, compared with Fig. 3, the temperature of the Fourier’ s slab increases faster while the temperature of the Cattaneo’s slab increases slower than that of the perfect thermal contact case. This is because the heat transfer from the Fourier’s region to the Cattaneo’s region is resisted by the imperfect interface in such a way more energy in the Fourier’s region is reserved to increase the temperature of the Slab 1 instead of transferring to the Slab 2. Besides, the constraints are excellent preserved in the machine error (1E-16). 3.3. Numerical example III : F 1T ¼ 0:1; F 2T ¼ 0:5 In this case, we consider the thermal contact of a general case emanating from the C-F heat conduction model. The heat conduc-
tivities are also considered to be different of the two slabs as: k1 ¼ 1 and k2 ¼ 2. Different heat conduction processes can be observed in Fig. 5. In general, the numerical results reveal the underlying physics of different F T . It’s worth noting that the parameter F T represents the portion of the Fourier-type in such a way the feature of telegraphic propagation in this case is stronger in the Slab-1 with smaller F T . In addition, Fig. 5(b) in the perfect thermal contact case shows that the constraints are automatically achieved at the beginning of the simulation. This may be due to that the heat signal is very weak at the interface such that the constraint does not play a role in the algorithm. Also, most of the heat is stored in the Slab-1 due to the thermal resistance, which leads the temperature of Slab-1 increases drastically while the heat of Slab-2 illustrates a gentle transfer. 4. Concluding remarks In this work, we have presented the formalism of a generalized heat conduction model with the consideration of perfect/imperfect thermal interface, which unifies both the Cattaneo’s and Jeffreys’ heat conduction (second-order time dependent system in terms of the single-field form) and the Fourier ’s type heat conduction (first order time dependent system in terms of the single-field form) into a simple two-field first-order time dependent differential algebraic framework in which the interface conditions are treated as algebraic equations, classified as Index 2 Constraints. The GS4 DAE Index 2 time integration is applied to carry out the numerical simulations. The numerical results from the GS4 DAE Index 2 time integration show the particular features of the thermal contacts between either Fourier- and/or Cattaneo- and/or Jeffreys- types with perfect/imperfect thermal interfaces. The proposed method has be proved to have the ability of describing
896
T. Xue et al. / International Journal of Heat and Mass Transfer 116 (2018) 889–896
the real heat conduction existing in the engineering problems such as multi-layer composites, film-substrate structures, and bio-tissue system, etc. A direct application of the proposed computational framework can be taken to is that the heat transport in the skins, which are usually comprised by epidermis, dermis, and hypodemis [20,21]. Each layer has its own material properties and relaxation times, which is a practical problem of the numerical example III. To model the heat process in such complex problem, the proposed computational framework can be directly implemented. Moreover, the structure of a general optical device is usually comprised by thin optical films and a thick subtract. The heat conduction of the film-subtract system has been described by the heat conduction model with relaxation time in the previous studies [9]. As a matter of fact, the thickness of the thin film may be at a close order of magnitude of its heat carrier’s mean free path, whereas the size of the subtract is usually much larger than the mean free path of its heat carriers, especially for thick metal subtracts, such as Fe and Cu. Using the same-type of heat conduction model seems to lack the reliability. The coupling different types of heat conduction models, such as the thermal contact between Fourier’ s and non-Fourier’ s heat conduction, appears to be a computational challenge in the previous studies since the heat flux continuity is difficult to be achieved. However, it has been proved to be very easy for the proposed computational framework and the numerical examples II and III have revealed this particular feature of the proposed method. The further works will mainly focus on the following two aspects: (i) the thermal contact in highly heterogeneous systems and the cases which are small in both the space and time; (ii) the application of general interface conditions on various imperfect thermal contacts. Acknowledgments We wish to thank the reviewers for helpful suggestions. We also would like to thank Prof. Kumar K. Tamma, Dr. Masao Shimada, and Allon Zhu from the University of Minnesota for their advice and support. Conflict of interest Authors declare that there is no conflict of interest. References [1] D. Hasselman, L.F. Johnson, Effective thermal conductivity of composites with interfacial thermal barrier resistance, J. Compos. Mater. 21 (6) (1987) 508–515.
[2] S. Torquato, M. Rintoul, Effect of the interface on the properties of composite media, Phys. Rev. Lett. 75 (22) (1995) 4067. [3] J. Yvonnet, Q.-C. He, Q.-Z. Zhu, J.-F. Shao, A general and efficient computational procedure for modelling the kapitza thermal resistance based on xfem, Comput. Mater. Sci. 50 (4) (2011) 1220–1224. [4] A. Javili, S. Kaessmair, P. Steinmann, General imperfect interfaces, Comput. Methods Appl. Mech. Eng. 275 (2014) 76–97. [5] A. Jain, K. Tamma, Parabolic heat conduction specialized applications involving imperfect contact surfaces: local discontinuous galerkin finite element method-part 2, J. Therm. Stresses 33 (4) (2010) 344–355. [6] H.-T. Chen, K.-C. Liu, Study of hyperbolic heat conduction problem in the film and substrate composite with the interface resistance, Jpn. J. Appl. Phys. 41 (10R) (2002) 6267. [7] D.Y. Tzou, Reflection and refraction of thermal waves from a surface or an interface between dissimilar materials, Int. J. Heat Mass Transfer 36 (2) (1993) 401–410. [8] K.-C. Liu, Analysis of dual-phase-lag thermal behaviour in layered films with temperature-dependent interface thermal resistance, J. Phys. D: Appl. Phys. 38 (19) (2005) 3722. [9] K.-C. Liu, Analysis of thermal behavior in multi-layer metal thin films based on hyperbolic two-step model, Int. J. Heat Mass Transfer 50 (7) (2007) 1397– 1407. [10] 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. [11] K.-C. Liu, H.-T. Chen, Analysis for the dual-phase-lag bio-heat transfer during magnetic hyperthermia treatment, Int. J. Heat Mass Transfer 52 (5) (2009) 1185–1192. [12] G. Yeoh, X. Gu, V. Timchenko, S. Valenzuela, B. Cornell, High order accurate dual-phase-lag numerical model for microscopic heating in multiple domains, Int. Commun. Heat Mass Transfer 78 (2016) 21–28. [13] W. Little, The transport of heat between dissimilar solids at low temperatures, Can. J. Phys. 37 (3) (1959) 334–349. [14] K.K. Tamma, X. Zhou, Macroscale and microscale thermal transport and thermo-mechanical interactions: some noteworthy perspectives, J. Therm. Stresses 21 (3–4) (1998) 405–449. [15] X. Zhou, K. Tamma, A new unified theory underlying time dependent linear first-order systems: a prelude to algorithms by design, Int. J. Numer. Meth. Eng. 60 (10) (2004) 1699–1740. [16] K. Tamma, D. Sha, X. Zhou, Time discretized operators. Part 1: Towards the theoretical design of a new generation of a generalized family of unconditionally stable implicit and explicit representations of arbitrary order for computational dynamics, Comput. Meth. Appl. Mech. Eng. 192 (3) (2003) 257–290. [17] S. Masuri, M. Sellier, X. Zhou, K. Tamma, Design of order-preserving algorithms for transient first-order systems with controllable numerical dissipation, Int. J. Numer. Meth. Eng. 88 (13) (2011) 1411–1448. [18] K.E. Jansen, C.H. Whiting, G.M. Hulbert, A generalized-a method for integrating the filtered navier–stokes equations with a stabilized finite element method, Comput. Meth. Appl. Mech. Eng. 190 (3) (2000) 305–319. [19] T. Xue, K.K. Tamma, X. Zhang, A consistent moving particle system simulation method: applications to parabolic/hyperbolic heat conduction type problems, Int. J. Heat Mass Transfer 101 (2016) 365–372. [20] F. Xu, K. Seffen, T. Lu, Non-fourier analysis of skin biothermomechanics, Int. J. Heat Mass Transfer 51 (9) (2008) 2237–2259. [21] J.T. Whitton, J. Everall, The thickness of the epidermis, Br. J. Dermatol. 89 (5) (1973) 467–476.