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

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 𝑞_𝑖,𝑖𝑞−𝑓=0inΩ,_𝑖+𝑘_𝑖𝑗𝑇_,𝑗=0inΩ, (1) 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 𝜕Ω_D and Neumann boundary conditions 𝜕Ω_N (i.e., 𝜕Ω_D∪𝜕Ω_N=𝜕Ω and 𝜕Ω_D∩𝜕Ω_N=∅): 𝑇=𝑇on𝜕Ω_D,𝑛_𝑖𝑞_𝑖=𝑞on𝜕Ω_N, (2) 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: _Ω𝑣𝑞_𝑖,𝑖−𝑓dΩ=0onΩ,_Ω𝑤_𝑖𝑗𝑞_𝑗+𝑘_𝑗𝑚𝑇_,𝑚dΩ=0onΩ. (3)

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 Ω^𝑒1∩Ω^𝑒2=∅, and it may have hanging nodes and elements of various shapes. Now the weak form of the problem can be given by _Ω𝑒𝑣^𝑒𝑞^𝑒_𝑖,𝑖−𝑓^𝑒dΩ=0,_Ω𝑒𝑤^𝑒_𝑖𝑗𝑞^𝑒_𝑗+𝑘^𝑒_𝑗𝑚𝑇^𝑒_,𝑚dΩ=0. (4)(5)

Integrating by parts the terms associated with the divergence in (4) and gradient in (5) lead to: _Ω𝑒𝑣^𝑒𝑞^𝑒_𝑖,𝑖dΩ=_𝜕Ω𝑒𝑣^𝑒𝑛^𝑒_𝑖dΓ−_Ω𝑒𝑣^𝑒_,𝑖𝑞^𝑒_𝑖dΩ,_Ω𝑒𝑤^𝑒_𝑖𝑗𝑘^𝑒_𝑗𝑚𝑇^𝑒_,𝑚dΩ=_𝜕Ω𝑒𝑇^𝑒𝑘^𝑒_𝑖𝑚𝑤^𝑒_𝑚𝑗𝑛^𝑒_𝑗−dΓ_Ω𝑒𝑇^𝑒𝑘^𝑒_𝑗𝑚𝑤^𝑒_{𝑚𝑖,𝑖}dΩ, (6)(7) where 𝑛^𝑒_𝑖 is the outward unit normal to the boundary 𝜕Ω^𝑒.

Substituting (6) into (4) and substituting (7) into (5), the following form can be obtained: _Ω𝑒𝑣^𝑒_,𝑖𝑞^𝑒_𝑖+𝑣^𝑒𝑓^𝑒dΩ=_𝜕Ω𝑒𝑣^𝑒𝑛^𝑒_𝑖𝑞^𝑒_𝑖dΓ,_Ω𝑒𝑇^𝑒𝑘^𝑒_𝑖𝑚𝑤^𝑒_{𝑚𝑗,𝑗}−𝑤^𝑒_𝑖𝑗𝑞^𝑒_𝑗dΩ=_𝜕Ω𝑒𝑇^𝑒𝑛^𝑒_𝑗𝑘^𝑒_𝑗𝑚𝑤^𝑒_𝑚𝑖dΓ. (8)

In (8) we note that the values of 𝑞^𝑒_𝑖 and 𝑇^𝑒 are required in the boundary integration terms ∮_𝜕Ω𝑒𝑣^𝑒𝑛^𝑒_𝑖𝑞^𝑒_𝑖dΓ and ∮_𝜕Ω𝑒𝑇^𝑒𝑛^𝑒_𝑗(𝑘^𝑒_𝑗𝑚𝑤^𝑒_𝑚𝑖)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 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: _Ω𝑒𝑣^𝑒_,𝑖𝑞^𝑒_𝑖+𝑣^𝑒𝑓^𝑒dΩ=_𝜕Ω𝑒𝑣^𝑒𝑛^𝑒_𝑖𝑞^𝑒_𝑖dΓ,_Ω𝑒𝑇^𝑒𝑘^𝑒_𝑖𝑚𝑤^𝑒_{𝑚𝑗,𝑗}−𝑤^𝑒_𝑖𝑗𝑞^𝑒_𝑗dΩ=_𝜕Ω𝑒𝑇^𝑒𝑛^𝑒_𝑗𝑘^𝑒_𝑗𝑚𝑤^𝑒_𝑚𝑖dΓ. (9)

In the present research, the numerical fluxes at the interior boundary Γ^𝑒𝑏 of element 𝑒 neighbored with element 𝑒𝑏 are defined by 𝑇^𝑒=𝛼_T𝑇^𝑒+𝛽_T𝑇^𝑒𝑏onΓ^𝑒𝑏,𝑞^𝑒_𝑖=𝛼_q𝑞^𝑒_𝑖+𝛽_q𝑞_𝑖^𝑒𝑏onΓ^𝑒𝑏. (10)

From a consistency point of view, the real valued numerical flux coefficients 𝛼_T, 𝛽_T, 𝛼_q, and 𝛽_q should satisfy the following constraint: 𝛼_T+𝛽_T𝛼=1,_q+𝛽_q=1. (11)

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𝑇^𝑒+𝑇^𝑒𝑏onΓ^𝑒𝑏,𝑞^𝑒_𝑖=12𝑞^𝑒_𝑖+𝑞_𝑖^𝑒𝑏onΓ^𝑒𝑏. (12) Numerical fluxes (12) implys that 𝛼_T=𝛽_T=𝛼_q=𝛽_𝑞=1/2.

If part of the element boundary 𝜕Ω^ext∈(𝜕Ω_D∪𝜕Ω_N), then the numerical fluxes on the external boundary 𝜕Ω^ext are defined by 𝑇^𝑒=𝑇,if𝜕Ω^ext∈𝜕Ω_D,𝑇^𝑒,if𝜕Ω^ext∈𝜕Ω_N,𝑛^𝑒_𝑖𝑞^𝑒_𝑖=𝑞,if𝜕Ω^ext∈𝜕Ω_N,𝑛^𝑒_𝑖𝑞^𝑒_𝑖,if𝜕Ω^ext∈𝜕Ω_D. (13)

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 _Ω𝑒𝑣^𝑒_,𝑖𝑞^𝑒_𝑖+𝑣^𝑒𝑓^𝑒dΩ+^𝑁𝑒𝑏_{𝑒𝑏=1Γ𝑒𝑏}𝜏𝑣^𝑒𝑇^𝑒−𝑇^𝑒𝑏=dΓ_𝜕Ω𝑒𝑣^𝑒𝑛^𝑒_𝑖𝑞^𝑒_𝑖dΓ,_Ω𝑒𝑇^𝑒𝑘^𝑒_𝑖𝑚𝑤^𝑒_{𝑚𝑗,𝑗}−𝑤^𝑒_𝑖𝑗𝑞^𝑒_𝑗dΩ=_𝜕Ω𝑒𝑇^𝑒𝑛^𝑒_𝑗𝑘^𝑒_𝑗𝑚𝑤^𝑒_𝑚𝑖dΓ, (14)(15) where 𝑁^𝑒_𝑏 denotes the number of boundary segments in element 𝑒, and Γ^𝑒𝑏 denotes the 𝑒𝑏th boundary segment in element 𝑒. The stabilization parameter 𝜏 should be 𝑂(‖𝑘^𝑒_𝑖𝑗‖/ℎ)>0, 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: 𝑇^𝑒=𝑁^𝑒_1𝑖ℑ^𝑒_𝑖,𝑞^𝑒_𝑖=𝑁^𝑒_2𝑖𝑗ℚ^𝑒_𝑗, (16) where ℑ^𝑒_𝑖 and ℚ^𝑒_𝑗 denote listings of nodal temperatures and heat fluxes for element 𝑒 and 𝑁^𝑒_1𝑖 and 𝑁^𝑒_2𝑖𝑗 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 𝑣^𝑒=𝑁^𝑒_1𝑖,𝑤^𝑒_𝑖𝑗=𝑁^𝑒_2𝑗𝑖. (17)

Substituting (16) and (17) into (14) and (15), the following discontinuous Galerkin finite element equations can be obtained 𝐾^𝑒_3𝑖𝑗−𝐾^𝑒_1𝑖𝑗ℚ^𝑒_𝑗+^𝑁𝑒𝑏_𝑒𝑏=1𝐾^𝑒𝑏_3𝑖𝑗ℚ_𝑗^𝑒𝑏−𝐾^𝑒_2𝑖𝑗ℑ^𝑒_𝑗−^𝑁𝑒𝑏_𝑒𝑏=1𝐾^𝑒𝑏_2𝑖𝑗ℑ_𝑗^𝑒𝑏=𝐹^𝑒_𝑇𝑖,𝐾^𝑒_5𝑖𝑗ℚ^𝑒_𝑗+𝐾^𝑒_6𝑖𝑗−𝐾^𝑒_4𝑖𝑗ℑ^𝑒_𝑗+^𝑁𝑒𝑏_𝑒𝑏=1𝐾^𝑒𝑏_6𝑖𝑗ℑ_𝑗^𝑒𝑏=0, (18) where we have used the definitions described in the following: 𝐾^𝑒_1𝑖𝑗=_Ω𝑒𝑁^𝑒_1𝑖,𝑚𝑁^𝑒_2𝑚𝑖𝐾dΩ,_{𝑒2𝑖𝑗}=_𝜕Ω𝑒𝜏𝑁^𝑒_1𝑖𝑁^𝑒_1𝑗𝐾dΩ,_{𝑒𝑏2𝑖𝑗}=_Γ𝑒𝑏−𝜏𝑁^𝑒_1𝑖𝑁^𝑒𝑏_1𝑗𝐾dΓ,_{𝑒3𝑖𝑗}=_𝜕Ω𝑒𝛼_q𝑁^𝑒_1𝑖𝑛^𝑒_𝑚𝑁^𝑒_2𝑚𝑗𝐾dΓ,_{𝑒𝑏3𝑖𝑗}=−_Γ𝑒𝑏𝛽_q𝑁^𝑒_1𝑖𝑛^𝑒_𝑚𝑁^𝑒𝑏_2𝑚𝑗𝐹dΓ,_𝑒𝑇𝑖=_Ω𝑒𝑁^𝑒_1𝑖𝑓^𝑒𝐾dΩ,_{𝑒4𝑖𝑗}=_Ω𝑒𝑁^𝑒_1𝑖𝑘^𝑒_𝑚𝑛𝑁^𝑒_{2𝑛𝑚,𝑗}𝐾dΩ,_{𝑒5𝑖𝑗}=_Ω𝑒−𝑁^𝑒_2𝑖𝑚𝑁^𝑒_2𝑚𝑗𝐾dΩ,_{𝑒6𝑖𝑗}=_𝜕Ω𝑒𝛼_T𝑁^𝑒_1𝑖𝑘^𝑒_𝑗𝑛𝑁^𝑒_2𝑛𝑚𝑛^𝑒_𝑚𝐾dΓ,_{𝑒𝑏6𝑖𝑗}=_Γ𝑒𝑏𝛽_T𝑁^𝑒𝑏_1𝑖𝑘^𝑒_𝑗𝑛𝑁^𝑒_2𝑛𝑚𝑛^𝑒_𝑚dΓ. (19)

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 ℚ^𝑒_𝑗=𝐾^𝑒_5𝑗𝑖⁻¹𝐾^𝑒_4𝑖𝑗−𝐾^𝑒_6𝑖𝑗ℑ^𝑒_𝑗−^𝑁𝑒𝑏_𝑒𝑏=1𝐾^𝑒𝑏_6𝑖𝑗ℑ_𝑗^𝑒𝑏. (20)

Then 𝑞^𝑒_𝑖 can be eliminated from the equations by substituting (20) into (14) to get the discontinuous Galerkin finite element equation on element domain: 𝐾^𝑒_𝑖𝑗ℑ^𝑒_𝑗+^𝑁𝑒𝑏_𝑒𝑏=1𝐾^𝑒𝑏_𝑖𝑗ℑ_𝑗^𝑒𝑏+𝐾^𝑒𝑏𝑏_𝑖𝑗ℑ_𝑗^𝑒𝑏𝑏=𝐹^𝑒_𝑇𝑖, (21) where 𝐾^𝑒_𝑖𝑗=𝐾^𝑒_3𝑖𝑗−𝐾^𝑒_1𝑖𝑗𝐾^𝑒_5𝑖𝑗⁻¹𝐾^𝑒_4𝑖𝑗−𝐾^𝑒_6𝑖𝑗−𝐾^𝑒_2𝑖𝑗,𝐾^𝑒𝑏_𝑖𝑗=𝐾^𝑒𝑏_3𝑖𝑗𝐾^𝑒_5𝑖𝑗⁻¹𝐾^𝑒𝑏_4𝑖𝑗−𝐾^𝑒𝑏_6𝑖𝑗−𝐾^𝑒_3𝑖𝑗−𝐾^𝑒_1𝑖𝑗𝐾^𝑒_5𝑖𝑗⁻¹𝐾^𝑒𝑏_6𝑖𝑗−𝐾^𝑒𝑏_2𝑖𝑗,𝐾^𝑒𝑏𝑏_𝑖𝑗=−^{𝑁𝑏𝑒𝑏𝑏}_{𝑒𝑏𝑏=1}𝐾^𝑒𝑏_3𝑖𝑗𝐾^𝑒𝑏_5𝑖𝑗⁻¹𝐾^𝑒𝑏𝑏_6𝑖𝑗. (22) 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 𝐾_𝑖𝑗ℑ_𝑗=𝐹_𝑖. (23)

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 ∇⋅𝝈+𝐟=𝟎onΩ, (24) where the operator ∇ is defined as 𝜕∇=𝜕𝜕𝑥000𝜕𝜕𝑧0𝜕𝜕𝑦0𝜕𝜕𝑦0𝜕𝜕𝑧𝜕𝜕𝑥00𝜕𝜕𝑧𝜕𝜕𝑦0𝜕𝑥. (25) The stress vector 𝝈 corresponds to the set 𝝈=𝜎_𝑥𝜎_𝑦𝜎_𝑧𝜏_𝑦𝑧𝜏_𝑧𝑥𝜏_𝑥𝑦^𝑇 and 𝐟 denotes the body force vector 𝐟=𝑓_𝑥𝑓_𝑦𝑓_𝑧^𝑇.

Let the boundary 𝜕Ω be decomposed into two parts: essential boundary 𝜕Ω_EB and natural boundary 𝜕Ω_NB, where 𝜕Ω_EB is fixed and a surface force 𝐩 is applied on 𝜕Ω_NB. That is, 𝐮=𝟎on𝜕Ω_EB,𝝈⋅𝐧=𝐩on𝜕Ω_NB. (26)

The weak form of the problem on an element Ω^𝑒 is given by _Ω𝑒𝐰^𝑒∇⋅𝝈^𝑒+𝐟^𝑒dΩ=0∀Ω^𝑒∈Ω, (27) where 𝐰^e denotes the weighting function.

Integrating by parts the terms associated with the operator in (27) leads to _𝜕Ω𝑒𝐰^𝑒𝐧^𝑒_𝑢⋅𝝈^𝑒dΓ−_Ω𝑒∇⋅𝐰^𝑒𝜎^𝑒dΩ+_Ω𝑒𝐰^𝑒𝐟^𝑒dΩ=0, (28) where 𝐧^𝑒=[𝑛^𝑒_𝑥𝑛^𝑒_𝑦𝑛^𝑒_𝑧]^𝑇 is the outward unit normal to the boundary 𝜕Ω^𝑒 and 𝐧^𝑒_𝑢=𝑛^𝑒_𝑥000𝑛^𝑒_𝑧𝑛^𝑒_𝑦0𝑛^𝑒_𝑦0𝑛^𝑒_𝑧0𝑛^𝑒_𝑥00𝑛^𝑒_𝑧𝑛^𝑒_𝑦𝑛^𝑒_𝑥0. (29)

Substituting (28) into (27), the following form can be obtained: _𝜕Ω𝑒𝐰^𝑒𝐧^𝑒_𝑢⋅𝝈^𝑒dΓ−_Ω𝑒∇⋅𝐰^𝑒𝝈^𝑒dΩ+_Ω𝑒𝐰^𝑒𝐟^𝑒dΩ=0. (30)

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 ∫_{Γ𝑖𝑒𝑏}𝐰^𝑒𝜏_𝑢(𝐮^𝑒−𝐮^𝑒𝑏−𝑔_𝑁𝐧^𝑒)dΓ 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: _Ω𝑒∇⋅𝐰^𝑒𝝈^𝑒dΩ+^𝑁𝑒_{𝑖=1Γ𝑖𝑒𝑏}𝐰^𝑒𝜏_𝑢𝐮^𝑒−𝐮^𝑒𝑏−𝑔_𝑁𝐧^𝑒dΓ=_Ω𝑒𝐰^𝑒𝐟^𝑒dΩ+_𝜕Ω𝑒𝐰^𝑒𝐧^𝑒_𝑢⋅𝝈^𝑒dΓ, (31) 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 𝝈^𝑒=12𝝈^𝑒+𝝈^𝑒𝑏onΓ_𝑖^𝑒𝑏. (32)

If part of the element boundary 𝜕Ω^ext∈(𝜕Ω_EB∪𝜕Ω_NB), then the stress numerical fluxes on the external boundary 𝜕Ω^ext are defined by 𝐧^𝑒_𝑢⋅𝝈^𝑒=𝐩,if𝜕Ω^ext∈𝜕Ω_NB,𝐧^𝑒_𝑢⋅𝝈^𝑒,if𝜕Ω^ext∈𝜕Ω_EB. (33)

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 _Ω𝑒∇⋅𝐰^𝑒𝝈^𝑒dΩ=_Ω𝑒𝐰^𝑒𝐟^𝑒dΩ+_𝜕Ω𝑒𝐰^𝑒𝐧^𝑒_𝑢⋅𝝈^𝑒dΓ. (34)

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: 𝐮^𝑒=𝐍^𝑒_𝑢𝐔^𝑒, (35) 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 𝐰^𝑒=𝐍^𝑒_𝑢^𝑇. (36)

The strain field of the element 𝑒 is given by 𝜺^𝑒=∇^𝑇𝐮^𝑒=𝐁^𝑒𝐔^𝑒 (37) with 𝐁^𝑒=∇^𝑇𝐍^𝑒_𝑢. (38)

Due to the assumption of linear elastic deformation, the relationship between the element stress and strain will be linear: 𝝈^𝑒=𝐃^𝑒⋅𝜺^𝑒−𝜶^𝑒𝑇^𝑒−𝑇₀, (39) 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 𝐊^𝑒_𝑢𝐔^𝑒+^𝑁𝑒_𝑖=1𝐊^𝑒_𝜏𝑢+𝐊^𝑒_𝑛𝑓𝑢𝐔^𝑒+𝐊^𝑒𝑏_𝜏𝑢+𝐊^𝑒𝑏_𝑛𝑓𝑢𝐔^𝑒𝑏=𝐅^𝑒_𝛼+𝐅^𝑒_𝑓+𝐅^𝑒_𝑁𝑢+^𝑁𝑒_𝑖=1𝐅^𝑒_𝑛𝑓+𝐅^𝑒𝑏_𝑛𝑓+𝐅^𝑒_gap, (40) where we have used the definitions described in the following: 𝐊^𝑒_𝑢=_Ω𝑒𝐁^𝑒𝑇𝐃^𝑒𝐁^𝑒𝐊dΩ,_𝑒𝜏𝑢=_{Γ𝑖𝑒𝑏}𝜏_𝑢𝐍^𝑒_𝑢^𝑇𝐍^𝑒_𝑢𝐊dΓ,_{𝑒𝑏𝜏𝑢}=−_Γeb𝑖𝜏_𝑢𝐍^𝑒_𝑢^𝑇𝐍_𝑢^𝑒𝑏𝐊dΓ,_{𝑒𝑛𝑓𝑢}=_{Γ𝑖𝑒𝑏}12𝐍^𝑒_𝑢^𝑇𝐧^𝑒_𝑢⋅𝐃^𝑒𝐁^𝑒𝐊dΓ,_{𝑒𝑏𝑛𝑓𝑢}=_{Γ𝑖𝑒𝑏}12𝐍^𝑒_𝑢^𝑇𝐧^𝑒_𝑢⋅𝐃^𝑒𝑏𝐁^𝑒𝑏𝐅dΓ,_𝑒𝛼=_Ω𝑒𝐁^𝑒𝑇𝐃^𝑒𝜶^𝑒𝑇^𝑒−𝑇₀dΩ^𝑒,𝐅^𝑒_𝑓=_Ω𝑒𝐍^𝑒_𝑢^𝑇𝐟^𝑒dΩ^𝑒,𝐅^𝑒_𝑛𝑓=_{Γ𝑖𝑒𝑏}12𝐍^𝑒_𝑢^𝑇𝐧^𝑒_𝑢⋅𝐃^𝑒𝜶^𝑒𝑇^𝑒−𝑇₀𝐅dΓ,_{𝑒𝑏𝑛𝑓}=_{Γ𝑖𝑒𝑏}12𝐍^𝑒_𝑢^𝑇𝐧^𝑒_𝑢⋅𝐃^𝑒𝑏𝜶^𝑒𝑏𝑇^𝑒𝑏−𝑇₀𝐅dΓ,_𝑒𝑁𝑢=_{Γ𝑖𝑒𝑏}𝐍^𝑒_𝑢^𝑇𝐅𝐩dΓ,_𝑒gap=_{Γ𝑖𝑒𝑏}𝜏_𝑢𝑔_𝑁𝐍^𝑒_𝑢^𝑇𝐧^𝑒dΓ. (41)

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)

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.

Figure 1 

Flowchart of computational algorithm in solving the thermomechanical coupling problem.

The convergence criterion is given by 𝛿_𝑇com=𝐓^(𝑛)−𝐓^(𝑛−1)𝐓^(𝑛)≤𝛿_𝑇,𝛿_𝑢com=𝐔^(𝑛)−𝐔^(𝑛−1)𝐔^(𝑛)≤𝛿_𝑢, (43) where 𝛿_𝑇com and 𝛿_𝑢com 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 𝐓^(𝑛−1)+𝛽⋅Δ𝐓^(𝑛) to update the current temperature vector instead of 𝐓^(𝑛−1)+Δ𝐓^(𝑛), 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. 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²), some fins (CDEF) are placed in contact with the plate. Sides EF and GF are convective boundaries with convection coefficient 3.6 W/m²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] 1𝑅=1.25𝑘^∗𝑚𝜎𝑝_𝑁𝑐₁1.617710⁶𝜎𝑚^{−𝑐20.95/(1+0.0711𝑐2)}, (44) where 𝑘^∗=2𝑘₁𝑘₂/(𝑘₁+𝑘₂) 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].

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 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)

Figure 7 

Temperature field of the structure with different interface pressures.

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.

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 cylinder 𝑟=0.1 m, wall thickness ℎ₁=0.04 m and ℎ₂=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 𝑇₂=300 K and the specified heat flux of the inner boundary of cylinder 1 is 𝑞₁=150 kWm². 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⁻⁸. 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)

Figure 9 

Numerical results of the structure with no initial gap: (a) temperature field; (b) radial displacement field.

(a)

(b)

Figure 10 

Stress fields of the structure with no initial gap: (a) radial stress; (b) hoop stress.

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.

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

Figure 12 

Computational model of heat-pipe-cooled leading edge thermal protection structure.

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.

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

(a)

(b)

Figure 14 

Structure response with no initial gap: (a) temperature field; (b) von Mises stress of the heat pipe.

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)

Figure 15 

The effect of thermal contact resistance with initial gaps: (a) stagnation temperature; (b) maximum von Mises stress.

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.

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.