# Applications of Discontinuous Galerkin Finite Element Method in Thermomechanical Coupling Problems with Imperfect Thermal Contact

Authors: Donghuan Liu; Yinghua Liu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Thermomechanical coupling problems with imperfect thermal contact are analyzed in the present paper with discontinuous Galerkin finite element method. The imperfect thermal contact condition is characterized by thermal contact resistance. The whole thermomechanical coupling problem is solved alternatively with the thermal subproblem and mechanical subproblem. Thermal contact resistance is introduced directly with the interface numerical flux of the present discontinuous Galerkin finite element method without using interface element as traditional continuous Galerkin finite element method does. Numerical results show the accuracy and feasibility of the present discontinuous Galerkin finite element method in solving thermomechanical coupling problems with imperfect thermal contact.


---

## Body

## 1. Introduction

Thermomechanical coupling problems exist in many engineering fields, such as hypersonic aircraft structures, nuclear fuel elements, heat exchanger, and metal forming. In these situations, predicting the temperature field as well as the stress field is of considerable applied importance. If different parts of the structure contact each other during the thermomechanical coupling process, then thermal contact resistance (TCR) appears, and such problems sometimes are called thermomechanical coupling contact problems.Other than the conventional thermomechanical coupling problems, thermal contact resistance plays an important role in the present analysis. Thermal contact resistance is caused by imperfect thermal contact; it arises in the contact region of two solids, as most of the heat through the interface flows through the actual contact spots, and this causes constrictions of heat flux lines through the bulk solid material of the contact interface, which leads to the thermal contact resistance. Over the preceding decades, many theoretical and experimental investigations have been performed to improve the knowledge of the thermal contact resistance. There have been several comprehensive reviews in thermal contact resistance modeling and experiment. Fieberg and Kneer developed an approach to derive the thermal contact resistance under high temperature and high pressure conditions based on transient infrared temperature measurements [1]. Shojaefard and Goudarzi proposed a numerical estimation of thermal contact resistance in contacting surfaces [2].  A predictive model for estimating thermal contact resistance between two nominal flat metallic rough surfaces had been developed and experimentally validated by Singhal et al. [3]. Liu and Shang proposed a physics-mechanism-based high temperature thermal contact conductance model [4]. Temizer and Wriggers developed a computational contact homogenization technique at the finite deformation regime to predict the macroscopic thermal response of contact interfaces between rough surface topographies [5]. Bahrami et al. reviewed the thermal contact resistance in a vacuum [6]; they divided the problem into three different parts: geometrical, mechanical, and thermal, and existing theories and models for each part are also reviewed. Because of the complexity of the thermal contact resistance, it is very hard to propose a unified model, so in the present research we will not focus on the thermal contact resistance model itself but the effect of thermal contact resistance on the thermomechanical coupling process.Thermal contact resistance directly leads to a temperature jump across the contact interface of the temperature filed, and then this changed temperature field will lead to the change of the stress field. As we know, thermal contact resistance is greatly dependent on contact stress, so the thermomechanical coupling problem occurred.The earliest work in this field can be traced back to Barber [8], who conducted a detailed investigation on the mechanism for thermoelastic instability (TEI). This work also laid a ground for later investigations on the phenomenon in both experiments and numerical simulations until recently [9–11]. It was found that the solution may become unstable when perfect contact condition is assumed at the contact interface [12], and in order to obtain a unique solution, one method is to introduce the pressure dependent thermal contact resistance. There are two main solution strategies for such thermomechanical coupling problems. The one is the monolithic schemes where temperature and displacement are solved together as a coupled system. The other is to treat the thermomechanical coupling problems as totally or partially uncoupled, and this strategy is usually permissible in practical applications. In the present paper, linear elastic small deformation assumption is adopted, and only the contribution of the temperature and pressure dependent thermal contact resistance is considered, so we can employ strategies known as staggered schemes. Therefore, the thermomechanical coupling problem is divided into two subproblems: a purely thermal subproblem at frozen configuration along with the thermal boundary conditions and a purely mechanical subproblem at frozen temperature along with the mechanical boundary conditions.In solving such thermomechanical coupling problems, the most important techniques are the implementation of thermal contact resistance and the impenetrability condition within the area of contact. Since the discontinuous Galerkin (DG) finite element method allows discontinuities of the physical unknowns within the interior of the problem domain, it seems to be a natural approach to capture the temperature jump at the interface in thermomechanical coupling contact problems. The discontinuous Galerkin finite element method results from the integration by parts on the finite element of the governing equations multiplied by discontinuous weighting functions. The first discontinuous Galerkin finite element method was introduced by Reed and Hill for the neutron transport equation [13], which is a linear hyperbolic partial differential equation. Although the original impetus of most researches on the discontinuous Galerkin method was the solution of hyperbolic and parabolic problems, the applications of the discontinuous Galerkin method have also spread to the elliptic problems [14]. Khalmanova and Costanzo proposed a space-time discontinuous Galerkin finite element method for fully coupled linear thermoelastodynamic problems with strain and heat flux discontinuities [15], but they did not consider the discontinuity caused by the imperfect contact on the material interface. Kanapady et al. provided the so-called local discontinuous Galerkin (LDG) method for solving various kinds of heat conduction problems with temperature discontinuity caused by thermal contact resistance [16, 17]. Liu et al. presented a discontinuous Galerkin finite element method for the heat conduction problems with local high gradient and thermal contact resistance [18]. They also proposed a discontinuous Galerkin finite element algorithm to solve the thermoelastic coupling problems caused by temperature and pressure dependent thermal contact resistance [19].The main purpose of the present paper is to propose a scheme of the discontinuous Galerkin finite element method and demonstrate its application in the modeling of thermomechanical coupling problems caused by imperfect thermal contact. The methods presented by Zienkiewicz et al. [20] and Zheng et al. [19] for one-dimensional and two-dimensional problems are also extended to three-dimensional situation.


## 2. Formulation of the Discontinuous Galerkin Finite Element Method

As mentioned before, thermomechanical coupling analysis presented here is divided into two subproblems, a purely thermal subproblem and a purely mechanical subproblem. Such an uncoupled formulation can be attractive for problems involving many degrees of freedom. This section firstly gives the discontinuous Galerkin finite element formulation of the thermal subproblem and mechanical subproblem and then shows the numerical implementation of the thermomechanical coupling analysis.


### 2.1. Thermal Analysis

Without loss of generality, we consider the following three-dimensional heat conduction problem on the domainΩ bounded by ∂Ω. The governing equations can be expressed as
(1)qi,i-f=0inΩ,qi+kijT,j=0inΩ,
where qi is the heat flux, f is the heat source, kij is thermal conductivity which may vary with temperature, and T is the temperature field. Here standard indicial notation is used, so repeated indices designate summations and a comma designates a partial derivative with respect to the indicated spatial variable.The boundary∂Ω is decomposed into a region of Dirichlet boundary conditions ∂ΩD and Neumann boundary conditions ∂ΩN (i.e., ∂ΩD∪∂ΩN=∂Ω and ∂ΩD∩∂ΩN=∅):
(2)T=T-on∂ΩD,niqi=q-on∂ΩN,
where T- and q- are the prescribed temperature and heat fluxs respectively, and ni is the outward unit normal to the boundary ∂Ω.In order to obtain the approximate solution of the problem, by taking the inner product of (1) with weighting functions v and wij, respectively, over the solution domain Ω we obtain the weak form of the problem:
(3)∫Ωv(qi,i-f)dΩ=0onΩ,∫Ωwij(qj+kjmT,m)dΩ=0onΩ.For a three-dimensional problem we assumeΞ(Ω) is a three-dimensional tessellation of the domain Ω. Let Ωe∈Ξ(Ω) bounded by ∂Ωe be a nonoverlapping element within the tessellation such that if e1≠e2 and then Ωe1∩Ωe2=∅, and it may have hanging nodes and elements of various shapes. Now the weak form of the problem can be given by
(4)∫Ωeve(qi,ie-fe)dΩ=0,(5)∫Ωewije(qje+kjmeT,me)dΩ=0.Integrating by parts the terms associated with the divergence in (4) and gradient in (5) lead to:
(6)∫Ωeveqi,iedΩ=∮∂ΩeveniedΓ-∫Ωev,ieqiedΩ,(7)∫Ωewije(kjmeT,me)dΩ=∮∂ΩeTe(kimewmje)njedΓ-∫ΩeTekjmewmi,iedΩ,
where nie is the outward unit normal to the boundary ∂Ωe.Substituting (6) into (4) and substituting (7) into (5), the following form can be obtained:
(8)∫Ωe(v,ieqie+vefe)dΩ=∮∂ΩevenieqiedΓ,∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeTenje(kjmewmie)dΓ.In (8) we note that the values of qie and Te are required in the boundary integration terms ∮∂ΩevenieqiedΓ and ∮∂ΩeTenje(kjmewmie)dΓ of each element, respectively. Due to the discontinuous nature of the approximation functions in the discontinuous Galerkin finite element, interelement discontinuity is allowed, so in the discontinuous Galerkin formulation continuity of the unknown variables is enforced by using stabilization terms and numerical fluxes between the interior elemental boundaries that will be discussed later.In the present discontinuous Galerkin formulation, we replace the element heat fluxqie and element temperature Te in the boundary integration terms by the so-called numerical fluxes q^ie and T^e, respectively, to obtain the discontinuous Galerkin finite element equations:
(9)∫Ωe(v,ieqie+vefe)dΩ=∮∂Ωevenieq^iedΓ,∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeT^enje(kjmewmie)dΓ.In the present research, the numerical fluxes at the interior boundaryΓeb of element e neighbored with element eb are defined by
(10)T^e=αTTe+βTTebonΓeb,q^ie=αqqie+βqqiebonΓeb.From a consistency point of view, the real valued numerical flux coefficientsαT, βT, αq, and βq should satisfy the following constraint:
(11)αT+βT=1,αq+βq=1.We can get different versions of discontinuous Galerkin finite element methodology by using different numerical flux coefficients. In the present research, the arithmetic average of the two values of temperature and heat flux at the boundary of the elements are employed as the following numerical fluxes:(12)T^e=12(Te+Teb)onΓeb,q^ie=12(qie+qieb)onΓeb.
Numerical fluxes (12) implys that αT=βT=αq=βq=1/2.If part of the element boundary∂Ωext∈(∂ΩD∪∂ΩN), then the numerical fluxes on the external boundary ∂Ωext are defined by
(13)T^e={T-,if∂Ωext∈∂ΩD,Te,if∂Ωext∈∂ΩN,nieq^ie={q-,if∂Ωext∈∂ΩN,nieqie,if∂Ωext∈∂ΩD.The stabilization terms in the formulation of discontinuous Galerkin method are of crucial importance. In the present research, we adopt the penalization of the jump of the temperature at the interface, so the stabilized discontinuous Galerkin finite element equations can be written as(14)∫Ωe(v,ieqie+vefe)dΩ+∑eb=1Nbe∫Γebτve(Te-Teb)dΓ=∮∂Ωevenieq^iedΓ,(15)∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeT^enje(kjmewmie)dΓ,
where Nbe denotes the number of boundary segments in element e, and  Γeb denotes the ebth boundary segment in element e. The stabilization parameter τ should be O(∥kije∥/h)>0, where h is an element size measure and ∥kije∥ is a norm of the thermal conduction coefficient matrix [20].Suppose that the temperatureTe and heat flux qie are approximated over a typical finite element Ωe by the following expressions:
(16)Te=N1ieℑie,qie=N2ijeℚje,
where ℑie and ℚje denote listings of nodal temperatures and heat fluxes for element e and N1ie and N2ije denote the expansion basis or shape functions. Here the shape functions can be constructed independently in each element domain Ωe, so that the order of the approximating polynomial can be easily changed from one element to the other and refinement of the grid can be achieved without taking into account the continuity restrictions which is typical of continuous finite element methods, and this makes the discontinuous Galerkin method  easily handle adaptivity strategies. In order to obtain a Galerkin style formulation, the weighting functions ve and wije are taken as
(17)ve=N1ie,wije=N2jie.Substituting (16) and (17) into (14) and (15), the following discontinuous Galerkin finite element equations can be obtained
(18)(K3ije-K1ije)ℚje+∑eb=1NbeK3ijebℚjeb-K2ijeℑje-∑eb=1NbeK2ijebℑjeb=FTie,K5ijeℚje+(K6ije-K4ije)ℑje+∑eb=1NbeK6ijebℑjeb=0,
where we have used the definitions described in the following:
(19)K1ije=∫ΩeN1i,meN2miedΩ,K2ije=∮∂ΩeτN1ieN1jedΩ,K2ijeb=∫Γeb-τN1ieN1jebdΓ,K3ije=∮∂ΩeαqN1ie(nmeN2mje)dΓ,K3ijeb=-∫ΓebβqN1ie(nmeN2mjeb)dΓ,FTie=∫ΩeN1iefedΩ,K4ije=∫ΩeN1iekmneN2nm,jedΩ,K5ije=∫Ωe-N2imeN2mjedΩ,K6ije=∮∂ΩeαTN1iekjneN2nmenmedΓ,K6ijeb=∫ΓebβTN1iebkjneN2nmenmedΓ.If the temperature numerical fluxT^e is independent of the heat flux on the element interior boundary, then the variable ℚje can be solved through (15) in an element-by-element fashion to get
(20)ℚje=(K5jie)-1((K4ije-K6ije)ℑje-∑eb=1Nbe(K6ijebℑjeb)).Thenqie can be eliminated from the equations by substituting (20) into (14) to get the discontinuous Galerkin finite element equation on element domain:
(21)Kijeℑje+∑eb=1Nbe(Kijebℑjeb+Kijebbℑjebb)=FTie,
where
(22)Kije=(K3ije-K1ije)(K5ije)-1(K4ije-K6ije)-K2ije,Kijeb=K3ijeb(K5ije)-1(K4ijeb-K6ijeb)-(K3ije-K1ije)(K5ije)-1K6ijeb-K2ijeb,Kijebb=-∑ebb=1Nbebb(K3ijeb(K5ijeb)-1K6ijebb).
Here the superscript “ebb” denotes the neighbor element of element eb.The assembly procedure for discontinuous Galerkin finite element method is the same as that used for continuous Galerkin finite element method, and the assembled system of equations corresponding to the unknown primary variables can be written as(23)Kijℑj=Fi.After the imposition of the boundary conditions procedure, the assembled equations can be solved by traditional techniques that are widely employed for the continuous Galerkin method; then we can use (15) to compute the heat flux or other desired quantities from the primary degrees of freedom computed from (23).


### 2.2. Mechanical Analysis

As in with the subsequent numerical implementation of the algorithm, in this subsection, the mechanical analysis process is given in matrix form. The governing equations in mechanical analysis can be expressed as(24)∇·σ+f=0onΩ,
where the operator ∇ is defined as
(25)∇=[∂∂x000∂∂z∂∂y0∂∂y0∂∂z0∂∂x00∂∂z∂∂y∂∂x0].
The stress vector σ corresponds to the set σ=[σxσyσzτyzτzxτxy]T  and f denotes the body force vector f=[fxfyfz]T.Let the boundary∂Ω be decomposed into two parts: essential boundary ∂ΩEB and natural boundary ∂ΩNB, where ∂ΩEB is fixed and a surface force p- is applied on ∂ΩNB. That is,
(26)u=0on∂ΩEB,σ·n=p-on∂ΩNB.The weak form of the problem on an elementΩe is given by
(27)∫Ωewe(∇·σe+fe)dΩ=0∀Ωe∈Ω,
where we denotes the weighting function.Integrating by parts the terms associated with the operator in (27) leads to
(28)∮∂Ωewe(nue·σe)dΓ-∫Ωe(∇·we)σedΩ+∫ΩewefedΩ=0,
where ne=[nxenyenze]T is the outward unit normal to the boundary ∂Ωe and
(29)nue=[nxe000nzenye0nye0nze0nxe00nzenyenxe0].Substituting (28) into (27), the following form can be obtained:
(30)∮∂Ωewe(nue·σe)dΓ-∫Ωe(∇·we)σedΩ+∫ΩewefedΩ=0.According to the classic treatment of discontinuous Galerkin finite element method [20], after replacing the boundary stress σe by the stress numerical flux σ^e and adding the stabilization terms ∫Γiebweτu(ue-ueb-gNne)dΓ to get weak enforcement of continuity, where gN denotes the gap function between element e and its neighbor eb, then the stabilized discontinuous Galerkin equation can be obtained:
(31)∫Ωe(∇·we)σedΩ+∑i=1Ne∫Γiebweτu(ue-ueb-gNne)dΓ=∫ΩewefedΩ+∮∂Ωewe(nue·σ^e)dΓ,
where Ne denotes the number of boundary segments in element e,Γieb denotes the ith common boundary segment of elements e and eb, and τu is the displacement stabilization parameter. Here the stress numerical flux σ^e is defined by
(32)σ^e=12(σe+σeb)onΓieb.If part of the element boundary∂Ωext∈(∂ΩEB∪∂ΩNB), then the stress numerical fluxes on the external boundary ∂Ωext are defined by
(33)nue·σ^e={p-,if∂Ωext∈∂ΩNB,nue·σe,if∂Ωext∈∂ΩEB.After the contact search process is finished, we just need to treat the continuity of the displacement of the contact segments, and here we adopt the penalty method to satisfy this geometrical constraint. That is to say, if contact occurs, then a stick mechanical contact condition needs to be established until the next iteration. In this case, the discontinuous Galerkin equation does not need to be modified; otherwise it implies that the gap is open and no mechanical contact occurs; then the stabilization terms must be eliminated from the discontinuous Galerkin equation; then we have(34)∫Ωe(∇·we)σedΩ=∫ΩewefedΩ+∮∂Ωewe(nue·σ^e)dΓ.It is worth noting that here we assume that once contact happens, it is in the so-called stick state, which implies that the nodes in contact are not allowed to move in both normal and tangential directions untill the next iteration to update the contact status again.Suppose that the element displacement fieldue is approximated by the following expression:
(35)ue=NueUe,
where Ue denotes listings of all nodal displacements for element e. Nue denotes the shape functions, and the shape functions can be constructed independently in each element domain.In order to obtain a Galerkin style formulation, the weighting functionswe are taken as
(36)we=(Nue)T.The strain field of the elemente is given by
(37)εe=∇Tue=BeUe
with
(38)Be=∇TNue.Due to the assumption of linear elastic deformation, the relationship between the element stress and strain will be linear:(39)σe=De·[εe-αe(Te-T0)],
where De is an elasticity matrix containing the appropriate material properties, T0 is the reference temperature, and αe is the thermal expansion coefficients matrix.Substituting (35)–(39) into (31), the following DG finite element equations can be obtained
(40)KueUe+∑i=1Ne((Kτue+Knfue)Ue+(Kτueb+Knfueb)Ueb)=Fαe+Ffe+FNue+∑i=1Ne(Fnfe+Fnfeb+Fgape),
where we have used the definitions described in the following:
(41)Kue=∫Ωe(Be)TDeBedΩ,Kτue=∫Γiebτu(Nue)TNuedΓ,Kτueb=-∫Γiebτu(Nue)TNuebdΓ,Knfue=∫Γieb12(Nue)Tnue·DeBedΓ,Knfueb=∫Γieb12(Nue)Tnue·DebBebdΓ,Fαe=∫Ωe(Be)TDeαe(Te-T0)dΩe,Ffe=∫Ωe(Nue)TfedΩe,Fnfe=∫Γieb12(Nue)Tnue·Deαe(Te-T0)dΓ,Fnfeb=∫Γieb12(Nue)Tnue·Debαeb(Teb-T0)dΓ,FNue=∫Γieb(Nue)Tp-dΓ,Fgape=∫ΓiebτugN(Nue)TnedΓ.For this problem, we also need to assemble the elemental discontinuous Galerkin equations to form the global equations, and after imposing the Dirichlet boundary conditions, the assembled equations can be solved by conventional techniques that are widely employed for the continuous Galerkin finite element method.


### 2.3. Numerical Implementation

Here we consider the continuous Galerkin finite element implementation as a baseline and comment on the essential modifications to obtain the discontinuous Galerkin one. Discontinuous Galerkin finite element method involves contributions of the stabilization and numerical flux integration terms on interior boundaries, which require the introduction of a loop over interior boundaries in addition to the loop over all the elements. This interaction between the two neighbored elements makes it necessary to implement a search algorithm that can detect the interface of the element and its neighbor element. In the present study, we give all the neighbor information such as neighbor elements and common sides when preparing the input gridding data.In the present research, the penalization of the jump of the temperature on the imperfect contact interface is eliminated to capture the temperature jump; at the same time the interface numerical flux of the heat flux is substituted by the definition of the thermal contact resistance:(42)ne·q^e=Teb-TeR.This treatment may lead to a novel approach to capture the temperature jump at the imperfect contact interface. Moreover, in the present paper the temperature and stress fields of each contact body are considered alternatively. In our study, since material properties are temperature dependent and thermal contact resistance depends on both of the temperature and contact pressure, the Picard iterative method is employed to solve the nonlinear equilibrium equations. The computational algorithm is given in Figure1.Figure 1
Flowchart of computational algorithm in solving the thermomechanical coupling problem.The convergence criterion is given by(43)δTcom=∥T(n)-T(n-1)∥∥T(n)∥≤δT,δucom=∥U(n)-U(n-1)∥∥U(n)∥≤δu,
where δTcom and δucom are the temperature and displacement convergence parameter computed, respectively, δT and δu are the temperature and displacement convergence tolerance, respectively, and the superscript “n” denotes the iteration number.In order to accelerate the convergence process and avoid the overmodification of the temperature vector, the relaxation technique is adopted in each iteration, which means to useT(n-1)+β·ΔT(n) to update the current temperature vector instead of T(n-1)+ΔT(n), where β is the relaxation parameter and 0<β≤1.


## 2.1. Thermal Analysis

Without loss of generality, we consider the following three-dimensional heat conduction problem on the domainΩ bounded by ∂Ω. The governing equations can be expressed as
(1)qi,i-f=0inΩ,qi+kijT,j=0inΩ,
where qi is the heat flux, f is the heat source, kij is thermal conductivity which may vary with temperature, and T is the temperature field. Here standard indicial notation is used, so repeated indices designate summations and a comma designates a partial derivative with respect to the indicated spatial variable.The boundary∂Ω is decomposed into a region of Dirichlet boundary conditions ∂ΩD and Neumann boundary conditions ∂ΩN (i.e., ∂ΩD∪∂ΩN=∂Ω and ∂ΩD∩∂ΩN=∅):
(2)T=T-on∂ΩD,niqi=q-on∂ΩN,
where T- and q- are the prescribed temperature and heat fluxs respectively, and ni is the outward unit normal to the boundary ∂Ω.In order to obtain the approximate solution of the problem, by taking the inner product of (1) with weighting functions v and wij, respectively, over the solution domain Ω we obtain the weak form of the problem:
(3)∫Ωv(qi,i-f)dΩ=0onΩ,∫Ωwij(qj+kjmT,m)dΩ=0onΩ.For a three-dimensional problem we assumeΞ(Ω) is a three-dimensional tessellation of the domain Ω. Let Ωe∈Ξ(Ω) bounded by ∂Ωe be a nonoverlapping element within the tessellation such that if e1≠e2 and then Ωe1∩Ωe2=∅, and it may have hanging nodes and elements of various shapes. Now the weak form of the problem can be given by
(4)∫Ωeve(qi,ie-fe)dΩ=0,(5)∫Ωewije(qje+kjmeT,me)dΩ=0.Integrating by parts the terms associated with the divergence in (4) and gradient in (5) lead to:
(6)∫Ωeveqi,iedΩ=∮∂ΩeveniedΓ-∫Ωev,ieqiedΩ,(7)∫Ωewije(kjmeT,me)dΩ=∮∂ΩeTe(kimewmje)njedΓ-∫ΩeTekjmewmi,iedΩ,
where nie is the outward unit normal to the boundary ∂Ωe.Substituting (6) into (4) and substituting (7) into (5), the following form can be obtained:
(8)∫Ωe(v,ieqie+vefe)dΩ=∮∂ΩevenieqiedΓ,∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeTenje(kjmewmie)dΓ.In (8) we note that the values of qie and Te are required in the boundary integration terms ∮∂ΩevenieqiedΓ and ∮∂ΩeTenje(kjmewmie)dΓ of each element, respectively. Due to the discontinuous nature of the approximation functions in the discontinuous Galerkin finite element, interelement discontinuity is allowed, so in the discontinuous Galerkin formulation continuity of the unknown variables is enforced by using stabilization terms and numerical fluxes between the interior elemental boundaries that will be discussed later.In the present discontinuous Galerkin formulation, we replace the element heat fluxqie and element temperature Te in the boundary integration terms by the so-called numerical fluxes q^ie and T^e, respectively, to obtain the discontinuous Galerkin finite element equations:
(9)∫Ωe(v,ieqie+vefe)dΩ=∮∂Ωevenieq^iedΓ,∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeT^enje(kjmewmie)dΓ.In the present research, the numerical fluxes at the interior boundaryΓeb of element e neighbored with element eb are defined by
(10)T^e=αTTe+βTTebonΓeb,q^ie=αqqie+βqqiebonΓeb.From a consistency point of view, the real valued numerical flux coefficientsαT, βT, αq, and βq should satisfy the following constraint:
(11)αT+βT=1,αq+βq=1.We can get different versions of discontinuous Galerkin finite element methodology by using different numerical flux coefficients. In the present research, the arithmetic average of the two values of temperature and heat flux at the boundary of the elements are employed as the following numerical fluxes:(12)T^e=12(Te+Teb)onΓeb,q^ie=12(qie+qieb)onΓeb.
Numerical fluxes (12) implys that αT=βT=αq=βq=1/2.If part of the element boundary∂Ωext∈(∂ΩD∪∂ΩN), then the numerical fluxes on the external boundary ∂Ωext are defined by
(13)T^e={T-,if∂Ωext∈∂ΩD,Te,if∂Ωext∈∂ΩN,nieq^ie={q-,if∂Ωext∈∂ΩN,nieqie,if∂Ωext∈∂ΩD.The stabilization terms in the formulation of discontinuous Galerkin method are of crucial importance. In the present research, we adopt the penalization of the jump of the temperature at the interface, so the stabilized discontinuous Galerkin finite element equations can be written as(14)∫Ωe(v,ieqie+vefe)dΩ+∑eb=1Nbe∫Γebτve(Te-Teb)dΓ=∮∂Ωevenieq^iedΓ,(15)∫Ωe(Tekimewmj,je-wijeqje)dΩ=∮∂ΩeT^enje(kjmewmie)dΓ,
where Nbe denotes the number of boundary segments in element e, and  Γeb denotes the ebth boundary segment in element e. The stabilization parameter τ should be O(∥kije∥/h)>0, where h is an element size measure and ∥kije∥ is a norm of the thermal conduction coefficient matrix [20].Suppose that the temperatureTe and heat flux qie are approximated over a typical finite element Ωe by the following expressions:
(16)Te=N1ieℑie,qie=N2ijeℚje,
where ℑie and ℚje denote listings of nodal temperatures and heat fluxes for element e and N1ie and N2ije denote the expansion basis or shape functions. Here the shape functions can be constructed independently in each element domain Ωe, so that the order of the approximating polynomial can be easily changed from one element to the other and refinement of the grid can be achieved without taking into account the continuity restrictions which is typical of continuous finite element methods, and this makes the discontinuous Galerkin method  easily handle adaptivity strategies. In order to obtain a Galerkin style formulation, the weighting functions ve and wije are taken as
(17)ve=N1ie,wije=N2jie.Substituting (16) and (17) into (14) and (15), the following discontinuous Galerkin finite element equations can be obtained
(18)(K3ije-K1ije)ℚje+∑eb=1NbeK3ijebℚjeb-K2ijeℑje-∑eb=1NbeK2ijebℑjeb=FTie,K5ijeℚje+(K6ije-K4ije)ℑje+∑eb=1NbeK6ijebℑjeb=0,
where we have used the definitions described in the following:
(19)K1ije=∫ΩeN1i,meN2miedΩ,K2ije=∮∂ΩeτN1ieN1jedΩ,K2ijeb=∫Γeb-τN1ieN1jebdΓ,K3ije=∮∂ΩeαqN1ie(nmeN2mje)dΓ,K3ijeb=-∫ΓebβqN1ie(nmeN2mjeb)dΓ,FTie=∫ΩeN1iefedΩ,K4ije=∫ΩeN1iekmneN2nm,jedΩ,K5ije=∫Ωe-N2imeN2mjedΩ,K6ije=∮∂ΩeαTN1iekjneN2nmenmedΓ,K6ijeb=∫ΓebβTN1iebkjneN2nmenmedΓ.If the temperature numerical fluxT^e is independent of the heat flux on the element interior boundary, then the variable ℚje can be solved through (15) in an element-by-element fashion to get
(20)ℚje=(K5jie)-1((K4ije-K6ije)ℑje-∑eb=1Nbe(K6ijebℑjeb)).Thenqie can be eliminated from the equations by substituting (20) into (14) to get the discontinuous Galerkin finite element equation on element domain:
(21)Kijeℑje+∑eb=1Nbe(Kijebℑjeb+Kijebbℑjebb)=FTie,
where
(22)Kije=(K3ije-K1ije)(K5ije)-1(K4ije-K6ije)-K2ije,Kijeb=K3ijeb(K5ije)-1(K4ijeb-K6ijeb)-(K3ije-K1ije)(K5ije)-1K6ijeb-K2ijeb,Kijebb=-∑ebb=1Nbebb(K3ijeb(K5ijeb)-1K6ijebb).
Here the superscript “ebb” denotes the neighbor element of element eb.The assembly procedure for discontinuous Galerkin finite element method is the same as that used for continuous Galerkin finite element method, and the assembled system of equations corresponding to the unknown primary variables can be written as(23)Kijℑj=Fi.After the imposition of the boundary conditions procedure, the assembled equations can be solved by traditional techniques that are widely employed for the continuous Galerkin method; then we can use (15) to compute the heat flux or other desired quantities from the primary degrees of freedom computed from (23).


## 2.2. Mechanical Analysis

As in with the subsequent numerical implementation of the algorithm, in this subsection, the mechanical analysis process is given in matrix form. The governing equations in mechanical analysis can be expressed as(24)∇·σ+f=0onΩ,
where the operator ∇ is defined as
(25)∇=[∂∂x000∂∂z∂∂y0∂∂y0∂∂z0∂∂x00∂∂z∂∂y∂∂x0].
The stress vector σ corresponds to the set σ=[σxσyσzτyzτzxτxy]T  and f denotes the body force vector f=[fxfyfz]T.Let the boundary∂Ω be decomposed into two parts: essential boundary ∂ΩEB and natural boundary ∂ΩNB, where ∂ΩEB is fixed and a surface force p- is applied on ∂ΩNB. That is,
(26)u=0on∂ΩEB,σ·n=p-on∂ΩNB.The weak form of the problem on an elementΩe is given by
(27)∫Ωewe(∇·σe+fe)dΩ=0∀Ωe∈Ω,
where we denotes the weighting function.Integrating by parts the terms associated with the operator in (27) leads to
(28)∮∂Ωewe(nue·σe)dΓ-∫Ωe(∇·we)σedΩ+∫ΩewefedΩ=0,
where ne=[nxenyenze]T is the outward unit normal to the boundary ∂Ωe and
(29)nue=[nxe000nzenye0nye0nze0nxe00nzenyenxe0].Substituting (28) into (27), the following form can be obtained:
(30)∮∂Ωewe(nue·σe)dΓ-∫Ωe(∇·we)σedΩ+∫ΩewefedΩ=0.According to the classic treatment of discontinuous Galerkin finite element method [20], after replacing the boundary stress σe by the stress numerical flux σ^e and adding the stabilization terms ∫Γiebweτu(ue-ueb-gNne)dΓ to get weak enforcement of continuity, where gN denotes the gap function between element e and its neighbor eb, then the stabilized discontinuous Galerkin equation can be obtained:
(31)∫Ωe(∇·we)σedΩ+∑i=1Ne∫Γiebweτu(ue-ueb-gNne)dΓ=∫ΩewefedΩ+∮∂Ωewe(nue·σ^e)dΓ,
where Ne denotes the number of boundary segments in element e,Γieb denotes the ith common boundary segment of elements e and eb, and τu is the displacement stabilization parameter. Here the stress numerical flux σ^e is defined by
(32)σ^e=12(σe+σeb)onΓieb.If part of the element boundary∂Ωext∈(∂ΩEB∪∂ΩNB), then the stress numerical fluxes on the external boundary ∂Ωext are defined by
(33)nue·σ^e={p-,if∂Ωext∈∂ΩNB,nue·σe,if∂Ωext∈∂ΩEB.After the contact search process is finished, we just need to treat the continuity of the displacement of the contact segments, and here we adopt the penalty method to satisfy this geometrical constraint. That is to say, if contact occurs, then a stick mechanical contact condition needs to be established until the next iteration. In this case, the discontinuous Galerkin equation does not need to be modified; otherwise it implies that the gap is open and no mechanical contact occurs; then the stabilization terms must be eliminated from the discontinuous Galerkin equation; then we have(34)∫Ωe(∇·we)σedΩ=∫ΩewefedΩ+∮∂Ωewe(nue·σ^e)dΓ.It is worth noting that here we assume that once contact happens, it is in the so-called stick state, which implies that the nodes in contact are not allowed to move in both normal and tangential directions untill the next iteration to update the contact status again.Suppose that the element displacement fieldue is approximated by the following expression:
(35)ue=NueUe,
where Ue denotes listings of all nodal displacements for element e. Nue denotes the shape functions, and the shape functions can be constructed independently in each element domain.In order to obtain a Galerkin style formulation, the weighting functionswe are taken as
(36)we=(Nue)T.The strain field of the elemente is given by
(37)εe=∇Tue=BeUe
with
(38)Be=∇TNue.Due to the assumption of linear elastic deformation, the relationship between the element stress and strain will be linear:(39)σe=De·[εe-αe(Te-T0)],
where De is an elasticity matrix containing the appropriate material properties, T0 is the reference temperature, and αe is the thermal expansion coefficients matrix.Substituting (35)–(39) into (31), the following DG finite element equations can be obtained
(40)KueUe+∑i=1Ne((Kτue+Knfue)Ue+(Kτueb+Knfueb)Ueb)=Fαe+Ffe+FNue+∑i=1Ne(Fnfe+Fnfeb+Fgape),
where we have used the definitions described in the following:
(41)Kue=∫Ωe(Be)TDeBedΩ,Kτue=∫Γiebτu(Nue)TNuedΓ,Kτueb=-∫Γiebτu(Nue)TNuebdΓ,Knfue=∫Γieb12(Nue)Tnue·DeBedΓ,Knfueb=∫Γieb12(Nue)Tnue·DebBebdΓ,Fαe=∫Ωe(Be)TDeαe(Te-T0)dΩe,Ffe=∫Ωe(Nue)TfedΩe,Fnfe=∫Γieb12(Nue)Tnue·Deαe(Te-T0)dΓ,Fnfeb=∫Γieb12(Nue)Tnue·Debαeb(Teb-T0)dΓ,FNue=∫Γieb(Nue)Tp-dΓ,Fgape=∫ΓiebτugN(Nue)TnedΓ.For this problem, we also need to assemble the elemental discontinuous Galerkin equations to form the global equations, and after imposing the Dirichlet boundary conditions, the assembled equations can be solved by conventional techniques that are widely employed for the continuous Galerkin finite element method.


## 2.3. Numerical Implementation

Here we consider the continuous Galerkin finite element implementation as a baseline and comment on the essential modifications to obtain the discontinuous Galerkin one. Discontinuous Galerkin finite element method involves contributions of the stabilization and numerical flux integration terms on interior boundaries, which require the introduction of a loop over interior boundaries in addition to the loop over all the elements. This interaction between the two neighbored elements makes it necessary to implement a search algorithm that can detect the interface of the element and its neighbor element. In the present study, we give all the neighbor information such as neighbor elements and common sides when preparing the input gridding data.In the present research, the penalization of the jump of the temperature on the imperfect contact interface is eliminated to capture the temperature jump; at the same time the interface numerical flux of the heat flux is substituted by the definition of the thermal contact resistance:(42)ne·q^e=Teb-TeR.This treatment may lead to a novel approach to capture the temperature jump at the imperfect contact interface. Moreover, in the present paper the temperature and stress fields of each contact body are considered alternatively. In our study, since material properties are temperature dependent and thermal contact resistance depends on both of the temperature and contact pressure, the Picard iterative method is employed to solve the nonlinear equilibrium equations. The computational algorithm is given in Figure1.Figure 1
Flowchart of computational algorithm in solving the thermomechanical coupling problem.The convergence criterion is given by(43)δTcom=∥T(n)-T(n-1)∥∥T(n)∥≤δT,δucom=∥U(n)-U(n-1)∥∥U(n)∥≤δu,
where δTcom and δucom are the temperature and displacement convergence parameter computed, respectively, δT and δu are the temperature and displacement convergence tolerance, respectively, and the superscript “n” denotes the iteration number.In order to accelerate the convergence process and avoid the overmodification of the temperature vector, the relaxation technique is adopted in each iteration, which means to useT(n-1)+β·ΔT(n) to update the current temperature vector instead of T(n-1)+ΔT(n), where β is the relaxation parameter and 0<β≤1.


## 3. Numerical Examples

### 3.1. Fin and Plate Problem

Considering the fin and plate problem that is typical in heat removal from electric devices, the goal of such structures is to keep the temperature of the plate as low as possible, and the key point is to keep the thermal contact resistance between the fin and the plate as low as possible. Figure2 shows the discretization of the problem. We consider a segment containing only a fin, and due to the symmetry of the problem only one-half of the geometry needs to be analyzed. To increase heat dissipation on the upper side (BG) of a plate (OABG) subjected to heat flux from the bottom OA (with prescribed heat flux 1000 W/m2), some fins (CDEF) are placed in contact with the plate. Sides EF and GF are convective boundaries with convection coefficient 3.6 W/m2K and convective exchange air at 293 K. Here we neglect the gas contribution and radiation effects to the thermal contact resistance and the TCR can be determined by [7]
(44)1R=1.25k*m-σ(pNc1(1.6177106σm-)-c2)0.95/(1+0.0711c2),
where k*=2k1k2/(k1+k2) is the mean thermal conductivity of the two contact materials, m- is the mean absolute asperity slope, σ is RMS surface roughness, pN is the apparent mechanical pressure, and c1, c2 are experimental hardness parameters determined with microhardness tests. The main interesting input material and structure data of the problem can be found in [7].Figure 2
Fin and plate problem description.Since the mechanical pressure between the fin and the plate is very important to the thermal contact resistance, a series of numerical tests have been carried out varying the mechanical pressure applied in the contact zone. The temperatures in some points of the structures, obtained by the present discontinuous Galerkin method, are compared with Wriggers and Zavarise’s results [7] shown in Figures 3, 4, 5, and 6.Figure 3
Comparisons of the variation of temperature at point A with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 4
Comparisons of the variation of temperature at point B with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 5
Comparisons of the variation of temperature at point C with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 6
Comparisons of the variation of temperature at point D with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].It can be observed from Figures3–6 that the present discontinuous Galerkin results are in general close to those obtained by Wriggers and Zavarise by using a two-node contact element [7], which shows the validation of the present method in modeling two-dimensional heat conduction problems with thermal contact resistance. Decreases of the temperature jump at the contact zone mean that heat can be easily transferred from the plate to the fin.Since high mechanical pressure may bring the structure failure at the contact zone, it is very important in engineering practices to determine the balance point of the mechanical pressure with the considerations both of the structure safety and the thermal dissipation efficiency, and the present discontinuous Galerkin method provides a novel approach to find this balance point efficiently and accurately. Figure7 shows the temperature field of the structure under different interface pressures.Temperature field of the structure with different interface pressures.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.2. Imperfect Thermal Contact of Two Cylinders

Here we consider two concentric cylinders in imperfect contact with each other, as shown in Figure8. Due to axial symmetry of the structure model, only a quarter of the whole model is considered in the numerical analysis.Figure 8
Illustration of the two concentric cylinders.Both cylinders are of identical materials and stress-free at 300 K. Plane stain conditions are assumed and the following data are employed in the analysis: inside radius of the inner cylinderr=0.1 m, wall thickness h1=0.04 m and h2=0.06 m (subscripts 1 and 2 refer to the inner and outer cylinders, resp.). The specified temperature of the outer boundary of cylinder 2 is T2=300 K and the specified heat flux of the inner boundary of cylinder 1 is q1=150 kWm2. The properties of the material and the thermal contact resistance model between the two cylinders and the input parameters of the TCR model can be found in [18].The computational algorithm described in the previous section is used to solve the present thermomechanical coupling problem, and both the temperature and displacement convergence tolerance are set to be 10−8. The computational results of temperature, radial stress, and hoop stress distributions in no initial gap case are illustrated in Figures 9-10, respectively.Numerical results of the structure with no initial gap: (a) temperature field; (b) radial displacement field.
(a)
(b)Stress fields of the structure with no initial gap: (a) radial stress; (b) hoop stress.
(a)
(b)As we can see from Figure10, since the interface contact pressure is over 200 MPa and the interface thermal contact resistance is very small, so the temperature jump at the interface of the two cylinders is very small. The present discontinuous Galerkin results are also compared well with the continuous Galerkin finite element method, as shown in Figure 11. Due to the different treatment of contact condition at the interface of the CG method and the present DG method, the discrepancy between two methods is larger around the middle part than both boundaries in as shown Figure 11.Figure 11
Radial stress distribution along radial direction for the two concentric cylinders.It is also noted that the radial and hoop stresses are very high; in order to reduce this thermal mismatch stress and smooth the stress distributions in the whole structure we usually set an initial radial gap between the two cylinders. This example shows the feasibility of the present method in solving two-dimensional thermomechanical coupling problems with thermal contact resistance.


### 3.3. Heat-Pipe-Cooled Leading Edge

The heat-pipe-cooled leading edge is an efficient thermal protection method in near space hypersonic aircraft design, and the heat transfer and thermomechanical coupling mechanisms are the key problems to be solved before its practical utility. Researches show that the difference between the experimental results and the theoretical results on the stagnation temperature of the heat-pipe-cooled leading edge is very large, and it is supposed that thermal contact resistance between heat pipe and high-temperature material around it has a significant influence on the thermal protection performances. The purpose of this numerical example is to investigate the effect of the thermal contact resistance to the thermomechanical coupling process of the heat-pipe-cooled leading edge with the present discontinuous Galerkin finite element method.The computational model of the heat-pipe-cooled leading edge thermal protection structure is illustrated in Figure12, and it includes a downwind area, an upwind area, a stagnation area, and a heat pipe inside it. Because of symmetry, half of the whole structure is considered here. The diameter of the heat pipe is 12.5 mm, and the thickness is 1.5 mm. The high temperature material around the heat pipe is three-dimensional braid C/C composite material with a thickness of 4.5 mm. The material of the heat pipe is Inconel 600, and the working fluid is alkali Na. The material properties can be found in [21].Figure 12
Computational model of heat-pipe-cooled leading edge thermal protection structure.The leading edge aerothermodynamic environment is given in Figure13. Radiation boundary condition is applied on the outer surface of the C/C material with a space temperature of 100 K. As heat pipe enables allowable equilibrium temperatures by providing exceptional effective thermal conductivity, we apply a forced boundary condition on the inner surface of the heat pipe with a constant temperature of 1000 K.Figure 13
Heat-pipe-cooled leading edge aerothermodynamic environment.As overestimate of the thermal contact resistance between the heat pipe and the C/C material can lead to a heavier thermal protection structure while underestimate of the thermal contact resistance can be a hidden danger to the structure, a rational thermal contact resistance model is the basis of the design and safety evaluation of the thermal protection structure. In the present paper, the thermal contact resistance formula is established based on the previous experimental research [21].Under the condition with no initial gap, the temperature field and von Mises stress field of the leading edge are given in Figure14.Structure response with no initial gap: (a) temperature field; (b) von Mises stress of the heat pipe.
(a)
(b)It can be seen from Figure14 that the stagnation point temperature is 1050 K, that is just 50 K higher than the temperature of the working fluid of the heat pipe and is totally under the capacity of the C/C material. However, the maximum von Mises stress of the heat pipe can up to 1150 MPa, which is far beyond the allowable temperature of the heat pipe material and is very dangerous. In order to reduce the thermal mismatch stress and smooth the stress distributions, an initial radial gap between the heat pipe and the C/C material is widely used.The effect of thermal contact resistance with different initial gaps is given in Figure15. Numerical results show that using heat-pipe-cooled leading edge can greatly reduce the maximum temperature to make sure that the stagnation temperature is in the tolerance zone (as shown in Figure 14), which is a very effective thermal protection concept. At the same time, using initial gap can greatly reduce the stress level of the structure, while increasing the interface thermal contact resistance and stagnation temperature (as shown in Figure 15(a)). With an initial gap of 65 μm, the maximum Mises stress of the structure is 62 MPa, but the maximum temperature is 1190 K. Structure strength and temperature should be considered simultaneously with the thermal contact resistance in structure safety evaluation of structures with initial gaps. If the initial gap is too small, the von Mises stress of the structure may cause structural failure because of the huge difference of thermal expansion coefficient of heat pipe and C/C material. If the initial gap is too big, the heat pipe and the C/C material may not contact each other and the intense heat of the stagnation area cannot redistribute quickly; then the whole structure may burn down. The feasibility zone of initial gap design is also given in Figure 16.The effect of thermal contact resistance with initial gaps: (a) stagnation temperature; (b) maximum von Mises stress.
(a)
(b)Figure 16
Feasibility zone of the initial gap design.


## 3.1. Fin and Plate Problem

Considering the fin and plate problem that is typical in heat removal from electric devices, the goal of such structures is to keep the temperature of the plate as low as possible, and the key point is to keep the thermal contact resistance between the fin and the plate as low as possible. Figure2 shows the discretization of the problem. We consider a segment containing only a fin, and due to the symmetry of the problem only one-half of the geometry needs to be analyzed. To increase heat dissipation on the upper side (BG) of a plate (OABG) subjected to heat flux from the bottom OA (with prescribed heat flux 1000 W/m2), some fins (CDEF) are placed in contact with the plate. Sides EF and GF are convective boundaries with convection coefficient 3.6 W/m2K and convective exchange air at 293 K. Here we neglect the gas contribution and radiation effects to the thermal contact resistance and the TCR can be determined by [7]
(44)1R=1.25k*m-σ(pNc1(1.6177106σm-)-c2)0.95/(1+0.0711c2),
where k*=2k1k2/(k1+k2) is the mean thermal conductivity of the two contact materials, m- is the mean absolute asperity slope, σ is RMS surface roughness, pN is the apparent mechanical pressure, and c1, c2 are experimental hardness parameters determined with microhardness tests. The main interesting input material and structure data of the problem can be found in [7].Figure 2
Fin and plate problem description.Since the mechanical pressure between the fin and the plate is very important to the thermal contact resistance, a series of numerical tests have been carried out varying the mechanical pressure applied in the contact zone. The temperatures in some points of the structures, obtained by the present discontinuous Galerkin method, are compared with Wriggers and Zavarise’s results [7] shown in Figures 3, 4, 5, and 6.Figure 3
Comparisons of the variation of temperature at point A with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 4
Comparisons of the variation of temperature at point B with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 5
Comparisons of the variation of temperature at point C with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].Figure 6
Comparisons of the variation of temperature at point D with the mechanical pressure at the contact zone between the present DG results and Wriggers and Zavarise’s results [7].It can be observed from Figures3–6 that the present discontinuous Galerkin results are in general close to those obtained by Wriggers and Zavarise by using a two-node contact element [7], which shows the validation of the present method in modeling two-dimensional heat conduction problems with thermal contact resistance. Decreases of the temperature jump at the contact zone mean that heat can be easily transferred from the plate to the fin.Since high mechanical pressure may bring the structure failure at the contact zone, it is very important in engineering practices to determine the balance point of the mechanical pressure with the considerations both of the structure safety and the thermal dissipation efficiency, and the present discontinuous Galerkin method provides a novel approach to find this balance point efficiently and accurately. Figure7 shows the temperature field of the structure under different interface pressures.Temperature field of the structure with different interface pressures.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2. Imperfect Thermal Contact of Two Cylinders

Here we consider two concentric cylinders in imperfect contact with each other, as shown in Figure8. Due to axial symmetry of the structure model, only a quarter of the whole model is considered in the numerical analysis.Figure 8
Illustration of the two concentric cylinders.Both cylinders are of identical materials and stress-free at 300 K. Plane stain conditions are assumed and the following data are employed in the analysis: inside radius of the inner cylinderr=0.1 m, wall thickness h1=0.04 m and h2=0.06 m (subscripts 1 and 2 refer to the inner and outer cylinders, resp.). The specified temperature of the outer boundary of cylinder 2 is T2=300 K and the specified heat flux of the inner boundary of cylinder 1 is q1=150 kWm2. The properties of the material and the thermal contact resistance model between the two cylinders and the input parameters of the TCR model can be found in [18].The computational algorithm described in the previous section is used to solve the present thermomechanical coupling problem, and both the temperature and displacement convergence tolerance are set to be 10−8. The computational results of temperature, radial stress, and hoop stress distributions in no initial gap case are illustrated in Figures 9-10, respectively.Numerical results of the structure with no initial gap: (a) temperature field; (b) radial displacement field.
(a)
(b)Stress fields of the structure with no initial gap: (a) radial stress; (b) hoop stress.
(a)
(b)As we can see from Figure10, since the interface contact pressure is over 200 MPa and the interface thermal contact resistance is very small, so the temperature jump at the interface of the two cylinders is very small. The present discontinuous Galerkin results are also compared well with the continuous Galerkin finite element method, as shown in Figure 11. Due to the different treatment of contact condition at the interface of the CG method and the present DG method, the discrepancy between two methods is larger around the middle part than both boundaries in as shown Figure 11.Figure 11
Radial stress distribution along radial direction for the two concentric cylinders.It is also noted that the radial and hoop stresses are very high; in order to reduce this thermal mismatch stress and smooth the stress distributions in the whole structure we usually set an initial radial gap between the two cylinders. This example shows the feasibility of the present method in solving two-dimensional thermomechanical coupling problems with thermal contact resistance.


## 3.3. Heat-Pipe-Cooled Leading Edge

The heat-pipe-cooled leading edge is an efficient thermal protection method in near space hypersonic aircraft design, and the heat transfer and thermomechanical coupling mechanisms are the key problems to be solved before its practical utility. Researches show that the difference between the experimental results and the theoretical results on the stagnation temperature of the heat-pipe-cooled leading edge is very large, and it is supposed that thermal contact resistance between heat pipe and high-temperature material around it has a significant influence on the thermal protection performances. The purpose of this numerical example is to investigate the effect of the thermal contact resistance to the thermomechanical coupling process of the heat-pipe-cooled leading edge with the present discontinuous Galerkin finite element method.The computational model of the heat-pipe-cooled leading edge thermal protection structure is illustrated in Figure12, and it includes a downwind area, an upwind area, a stagnation area, and a heat pipe inside it. Because of symmetry, half of the whole structure is considered here. The diameter of the heat pipe is 12.5 mm, and the thickness is 1.5 mm. The high temperature material around the heat pipe is three-dimensional braid C/C composite material with a thickness of 4.5 mm. The material of the heat pipe is Inconel 600, and the working fluid is alkali Na. The material properties can be found in [21].Figure 12
Computational model of heat-pipe-cooled leading edge thermal protection structure.The leading edge aerothermodynamic environment is given in Figure13. Radiation boundary condition is applied on the outer surface of the C/C material with a space temperature of 100 K. As heat pipe enables allowable equilibrium temperatures by providing exceptional effective thermal conductivity, we apply a forced boundary condition on the inner surface of the heat pipe with a constant temperature of 1000 K.Figure 13
Heat-pipe-cooled leading edge aerothermodynamic environment.As overestimate of the thermal contact resistance between the heat pipe and the C/C material can lead to a heavier thermal protection structure while underestimate of the thermal contact resistance can be a hidden danger to the structure, a rational thermal contact resistance model is the basis of the design and safety evaluation of the thermal protection structure. In the present paper, the thermal contact resistance formula is established based on the previous experimental research [21].Under the condition with no initial gap, the temperature field and von Mises stress field of the leading edge are given in Figure14.Structure response with no initial gap: (a) temperature field; (b) von Mises stress of the heat pipe.
(a)
(b)It can be seen from Figure14 that the stagnation point temperature is 1050 K, that is just 50 K higher than the temperature of the working fluid of the heat pipe and is totally under the capacity of the C/C material. However, the maximum von Mises stress of the heat pipe can up to 1150 MPa, which is far beyond the allowable temperature of the heat pipe material and is very dangerous. In order to reduce the thermal mismatch stress and smooth the stress distributions, an initial radial gap between the heat pipe and the C/C material is widely used.The effect of thermal contact resistance with different initial gaps is given in Figure15. Numerical results show that using heat-pipe-cooled leading edge can greatly reduce the maximum temperature to make sure that the stagnation temperature is in the tolerance zone (as shown in Figure 14), which is a very effective thermal protection concept. At the same time, using initial gap can greatly reduce the stress level of the structure, while increasing the interface thermal contact resistance and stagnation temperature (as shown in Figure 15(a)). With an initial gap of 65 μm, the maximum Mises stress of the structure is 62 MPa, but the maximum temperature is 1190 K. Structure strength and temperature should be considered simultaneously with the thermal contact resistance in structure safety evaluation of structures with initial gaps. If the initial gap is too small, the von Mises stress of the structure may cause structural failure because of the huge difference of thermal expansion coefficient of heat pipe and C/C material. If the initial gap is too big, the heat pipe and the C/C material may not contact each other and the intense heat of the stagnation area cannot redistribute quickly; then the whole structure may burn down. The feasibility zone of initial gap design is also given in Figure 16.The effect of thermal contact resistance with initial gaps: (a) stagnation temperature; (b) maximum von Mises stress.
(a)
(b)Figure 16
Feasibility zone of the initial gap design.


## 4. Conclusions

Thermomechanical coupling analysis with imperfect thermal contact is performed with the discontinuous Galerkin finite element method in the present paper. The coupling problem is divided into a thermal problem and a mechanical problem, and they are formulated with the discontinuous Galerkin finite element method, respectively. The imperfect thermal contact condition is realized directly at the contact boundary where thermal numerical flux is given by the definition of thermal contact resistance without using interface element in continuous Galerkin finite element method. Three numerical examples are given to show the feasibility and accuracy of the present discontinuous Galerkin finite element method in modeling with thermomechanical coupling problems with imperfect thermal contact.It also should be mentioned that there is still a lot of work to do for discontinuous Galerkin finite element method in further research. Firstly, the discontinuous nature of the shape functions introduces additional unknowns relative to the continuous Galerkin method, as we need to distinguish between values from different elements at the same node. Secondly, the implementation is more complicated than that of standard finite element procedures, especially in three-dimensional cases. This is evident in the treatment of the interfaces where information from neighbor elements is needed to construct the flux terms. However, the method promises of efficient parallel implementation and adaptivity, which can reduce the effect of the increase of degrees of freedom. Moreover, as the high gradient and discontinuities are always located in some parts of the mesh, the coupled discontinuous/continuous Galerkin method is more preferable [22]. So future work aims at further enhancing the efficiency of the computations, especially in regions where the solution is smooth.



---

*Source: 861417-2013-10-01.xml*