#### 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.

#### 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 where is the heat flux, is the heat source, is thermal conductivity which may vary with temperature, and 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 and Neumann boundary conditions (i.e., and ): where and are the prescribed temperature and heat fluxs respectively, and 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 and , respectively, over the solution domain we obtain the weak form of the problem:

For a three-dimensional problem we assume is a three-dimensional tessellation of the domain . Let bounded by be a nonoverlapping element within the tessellation such that if and then , and it may have hanging nodes and elements of various shapes. Now the weak form of the problem can be given by

Integrating by parts the terms associated with the divergence in (4) and gradient in (5) lead to: where is the outward unit normal to the boundary .

Substituting (6) into (4) and substituting (7) into (5), the following form can be obtained:

In (8) we note that the values of and are required in the boundary integration terms and 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 flux and element temperature in the boundary integration terms by the so-called numerical fluxes and , respectively, to obtain the discontinuous Galerkin finite element equations:

In the present research, the numerical fluxes at the interior boundary of element neighbored with element are defined by

From a consistency point of view, the real valued numerical flux coefficients , , , and should satisfy the following constraint:

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: Numerical fluxes (12) implys that .

If part of the element boundary , then the numerical fluxes on the external boundary are defined by

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 where denotes the number of boundary segments in element , and denotes the th boundary segment in element . The stabilization parameter should be , where is an element size measure and is a norm of the thermal conduction coefficient matrix [20].

Suppose that the temperature and heat flux are approximated over a typical finite element by the following expressions: where and denote listings of nodal temperatures and heat fluxes for element and and denote the expansion basis or shape functions. Here the shape functions can be constructed independently in each element domain , 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 and are taken as

Substituting (16) and (17) into (14) and (15), the following discontinuous Galerkin finite element equations can be obtained where we have used the definitions described in the following:

If the temperature numerical flux is independent of the heat flux on the element interior boundary, then the variable can be solved through (15) in an element-by-element fashion to get

Then can be eliminated from the equations by substituting (20) into (14) to get the discontinuous Galerkin finite element equation on element domain: where Here the superscript “” denotes the neighbor element of element .

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

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 where the operator is defined as The stress vector corresponds to the set and denotes the body force vector .

Let the boundary be decomposed into two parts: essential boundary and natural boundary , where is fixed and a surface force is applied on . That is,

The weak form of the problem on an element is given by where denotes the weighting function.

Integrating by parts the terms associated with the operator in (27) leads to where is the outward unit normal to the boundary and

Substituting (28) into (27), the following form can be obtained:

According to the classic treatment of discontinuous Galerkin finite element method [20], after replacing the boundary stress by the stress numerical flux and adding the stabilization terms to get weak enforcement of continuity, where denotes the gap function between element and its neighbor , then the stabilized discontinuous Galerkin equation can be obtained: where denotes the number of boundary segments in element , denotes the th common boundary segment of elements and , and is the displacement stabilization parameter. Here the stress numerical flux is defined by

If part of the element boundary , then the stress numerical fluxes on the external boundary are defined by

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

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 field is approximated by the following expression: where denotes listings of all nodal displacements for element . 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 functions are taken as

The strain field of the element is given by with

Due to the assumption of linear elastic deformation, the relationship between the element stress and strain will be linear: where is an elasticity matrix containing the appropriate material properties, is the reference temperature, and is the thermal expansion coefficients matrix.

Substituting (35)–(39) into (31), the following DG finite element equations can be obtained where we have used the definitions described in the following:

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:

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 Figure 1.

The convergence criterion is given by where and are the temperature and displacement convergence parameter computed, respectively, and are the temperature and displacement convergence tolerance, respectively, and the superscript “” 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 use to update the current temperature vector instead of , where is the relaxation parameter and .

#### 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. Figure 2 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/m^{2}), some fins (CDEF) are placed in contact with the plate. Sides EF and GF are convective boundaries with convection coefficient 3.6 W/m^{2}K 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]
where is the mean thermal conductivity of the two contact materials, is the mean absolute asperity slope, is RMS surface roughness, is the apparent mechanical pressure, and , are experimental hardness parameters determined with microhardness tests. The main interesting input material and structure data of the problem can be found in [7].

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.

It can be observed from Figures 3–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. Figure 7 shows the temperature field of the structure under 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 Figure 8. Due to axial symmetry of the structure model, only a quarter of the whole model is considered in the numerical analysis.

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 cylinder m, wall thickness m and m (subscripts 1 and 2 refer to the inner and outer cylinders, resp.). The specified temperature of the outer boundary of cylinder 2 is K and the specified heat flux of the inner boundary of cylinder 1 is kWm^{2}. 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.

**(a)**

**(b)**

**(a)**

**(b)**

As we can see from Figure 10, 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.

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 Figure 12, 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].

The leading edge aerothermodynamic environment is given in Figure 13. 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.

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 Figure 14.

**(a)**

**(b)**

It can be seen from Figure 14 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 Figure 15. 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.

**(a)**

**(b)**

#### 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.

#### Acknowledgments

This project is supported by the National High Technology Research and Development Program of China (Grant no. 2012AA03A513), the National Natural Science Foundation of China (Grant nos. 11302023 and 91016026), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant no. 20120006120008), and the Fundamental Research Funds for the Central Universities (Grant no. FRF-BR-13-003) which are greatly acknowledged.