Computational Mechanics https://doi.org/10.1007/s00466-018-1555-z

ORIGINAL PAPER

Application of an enriched FEM technique in thermo-mechanical contact problems A. R. Khoei1 · B. Bahmani1 Received: 28 July 2017 / Accepted: 30 January 2018 © Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract In this paper, an enriched FEM technique is employed for thermo-mechanical contact problem based on the extended finite element method. A fully coupled thermo-mechanical contact formulation is presented in the framework of X-FEM technique that takes into account the deformable continuum mechanics and the transient heat transfer analysis. The Coulomb frictional law is applied for the mechanical contact problem and a pressure dependent thermal contact model is employed through an explicit formulation in the weak form of X-FEM method. The equilibrium equations are discretized by the Newmark time splitting method and the final set of non-linear equations are solved based on the Newton–Raphson method using a staggered algorithm. Finally, in order to illustrate the capability of the proposed computational model several numerical examples are solved and the results are compared with those reported in literature. Keywords Enriched FEM · Thermo-mechanics · Frictional contact · Thermal contact

1 Introduction Thermo-mechanical modeling of contact problem is an essential tool in the design of heated bodies in order to achieve a proper resistance against thermal defects. The heat transfer analysis involving mechanical effects is interested in various engineering areas, such as metal forming process, crashworthiness and projectile impact, fault interaction in hot crust, and etc [1–4]. The presence of thermal cracks is highly feasible in such problems. Hence, a thermomechanical analysis of the cracked body plays a vital role in industrial instruments. One of the most important aspects in the thermo-mechanical analysis of a cracked domain is the study of thermal-contact behavior at the area of contact edges. The analytical solutions for contact problems are limited to special cases which inevitably lead to utilizing the numerical simulations. The FEM has been employed by many researchers to model the contact problem based on the minimization of a functional, which describes the governing equation of physical problem with special constraints,

B 1

A. R. Khoei [email protected] Department of Civil Engineering, Center of Excellence in Structures and Earthquake Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran

named as contact constraints. The contact constraints can be imposed into the weak formulation by various techniques, such as the penalty, Lagrange, augmented–Lagrange, mortar, Nitsche, and etc. However, the standard FEM method has some difficulties in modeling contact problems, such as the alignment of FE mesh with the contact surface and the establishment of physical connection between two contact bodies in numerical treatment. As a result, the extended finite element method (X-FEM) has been employed with the help of level set strategy to deal with these restrictions. The X-FEM was first developed by Belytschko and Black [5] to model the crack discontinuity and its propagation in the solid domain. The technique benefits from an appropriate space of approximation to capture individual characteristics of local physical phenomena, such as the singular behavior, sharp gradients, and strong discontinuities [6–8]. Implementation of the X-FEM in modeling of contact problem was performed by Dolbow et al. [9] using the LATIN method. Khoei and Nikbakht [10] presented the frictional contact behavior based on the penalty method within the X-FEM framework. The technique was then extended by Liu and Borja [11] to capture the path dependent crack propagation. Khoei and Vahab [12] proposed a numerical contact algorithm for simulation of saturated porous media with the X-FEM. Further research works were performed by researchers for better handling contact constrains in the

123

Computational Mechanics

X-FEM to eliminate the non-physical oscillations in numerical solution [13–16]. The X-FEM was first employed in thermo-mechanical problems by Duflot [17] to model the thermo-elastic fracture mechanics. Zamani and Eslami [18] developed an X-FEM model for crack propagation in thermoelasto-dynamic fracture problems. Khoei et al. [19] applied the X-FEM in thermo-hydro-mechanical modeling of impermeable discontinuity in saturated porous media. Shao et al. [20] proposed the X-FEM in thermo-mechanical analysis of cracked porous media considering the effects of fluid flow and heat transfer. Gill and Davey [21] employed an X-FEM thermo-mechanical analysis to model leaks through cracks using the enriched spaces for thermal and mechanical approximations. The thermal resistance at the contact interface between two bodies plays a crucial role in modeling the thermomechanical contact problem; so a proper constitutive law must be employed to represent the thermal behavior at contacting bodies. Yvonnet et al. [22] proposed the linear theory of heat transfer to model the thermal resistance by the XFEM, in which the effect of highly conductive interfaces was studied in the composite materials. Yvonnet et al. [23] modeled a special kind of thermal resistance interfaces, named Kapitza interface using the X-FEM. Gu et al. [24] developed a general imperfect interface model for thermal conduction in composites; the technique was then employed by Liu et al. [25] within the X-FEM. Javili et al. [26,27] studied the thermal resistance of different contact interfaces using the conventional thermo-mechanical FEM method. The thermal conduction in imperfect contact interfaces were studied by Jain and Tamma [28] and Gu et al. [29] for linear coupled multi-field phenomena. The thermo-elastic contact problem was first formulated within the standard FEM by Curnier and Taylor [30]. The coupled thermo-mechanical contact analysis was presented by Zavarise et al. [31] for real contact mechanisms. A thermo-mechanical contact friction model was developed by Wriggers and Miehe [32] and Zavarise et al. [33] in large deformation problems. A thermo-mechanical contact model considering plasticity behavior was employed by Pantuso et al. [34]. An adaptive FEM analysis was applied by Rieger and Wriggers [35] for thermo-mechanical coupled contact problems. A methodology to deal with the high temperature gradients was developed by Zavarise et al. [36] for frictional heating with application to penetration process. Hansen [37] proposed an approach based on the Mortar technique to impose the thermo-mechanical contact constraints using the Lagrange method. An advanced higher-order thermo-mechanical contact scheme that is based on a phase-field approach in a mortar contact framework was proposed by Hesch et al. [38]. Belghith et al. [39] developed a method based on the homogenization technique to model thermo-mechanical behavior of rough

123

contacting surfaces. Dittmann et al. [40] introduced a strategy for thermo-mechanical mortar contact problems based on the isogeometric analysis. Murashov and Panin [41] studied the effect of contact rough surfaces in thermomechanical contact analysis of work hardening contact interfaces. To the knowledge of authors, there is no explicit algorithm for modeling the thermo-mechanical contact problem based on an enriched FEM method. The aim of this research is to present a formulation for the pressure-dependent thermomechanical contact problem using the X-FEM method, in which the frictional sliding is modeled using the Coulomb friction law and the quasi-static infinitesimal deformation is adopted for problems presented here. In the thermomechanical contact model, the transient heat transfer analysis is performed, and the thermal contact condition is explicitly specified as a nonlinear constitutive law that describes the heat flux through the contact zone. The mechanical and thermal behaviors have two sources of coupling; firstly, a global source through the thermal expansion resulting from the temperature effect on mechanical deformation, and secondly which is a highly nonlinear local source through the pressure dependent thermal-contact model resulting from the contact forces that affect on the thermal-contact constitutive law. Finally, a staggered scheme is employed using the Newton–Raphson iterative procedure to solve the set of coupled nonlinear equations. The paper is organized as follows; in Sect. 2, the weak form of a fractured thermo-elastic domain is derived by introducing the internal boundary conditions on crack surfaces within the thermo-mechanical governing equations. In Sect. 3, the mechanical and thermal contact formulations are presented by imposing the contact constraints based on the penalty method. In Sect. 4, the system of nonlinear equations is described by employing the X-FEM technique for the spatial discretization and the Newmark scheme for the time domain discretization. In Sects. 5 and 6, the computational algorithm are presented for governing equations of the thermo-mechanical contact problem. Finally in Sect. 7, several numerical examples are solved to illustrate the capability of the proposed computational algorithm.

2 Governing equations of a fractured thermo-elastic medium In order to derive the governing equations of a cracked body in thermo-mechanical loading conditions, the weak form of governing equations is obtained within the X-FEM context employing the Divergence theorem. It must be noted that description of the deformation is restricted to infinitesimal strain and the linear elastic behavior is considered for the bulk material. Consider an arbitrary body  is bounded by

Computational Mechanics

Fig. 1 A cracked body in thermo-mechanical contact condition; problem definition

the external boundary , which is subjected to boundary conditions, as shown in Fig. 1. The mechanical boundary conditions are prescribed by the traction t on t and the prescribed displacement u˜ on u , such that t ∪ u =  and t ∩ u = ∅; the thermal boundary conditions are prescribed by the heat flux qn on q and the temperature θ˜ on θ with q ∪ θ =  and q ∩ θ = ∅. Moreover, the fractured domain contains an internal discontinuous interface denoted by d that generally includes the contact surface c . The well-known momentum balance equilibrium equation of a fractured body in quasi-static condition can be written as ∇ ·σ +ρb=0

(1)

where b is the body force per unit mass vector, ρ is the density of the material, and σ is the Cauchy stress tensor. The stress tensor is related to the strain tensor ε using an appropriate constitutive law. The above equilibrium equation is associated with the external and internal boundary conditions given by σ.n = t u = u˜ σ . nd =

tcont

on

t

on

u

on

c ⊂ d

(2)

where tcont is the contact traction imposed on the discontinuity surface d , n is the unit normal vector on the external boundary, and nd is the unit normal vector on the internal boundary from − to + . It must be noted that tcont has the non-zero value on the contact surface c and zero elsewhere because of the traction free boundary conditions on the opened discontinuities. The constitutive equation relating the stress tensor to the strain tensor is declared by the Duhamel-Neumann type of

Hooke’s law as σ = D : (ε − εθ ), where D is the elasticity matrix of the bulk material, and the strain tensor is defined as ε(u) = ∇ s u with ∇ s denoting the symmetrical spatial gradient operator, i.e. ∇ s = 21 (∇ +∇ T ). In the above definition, εθ is the thermal strain tensor defined by εθ (θ ) = α (θ − θ0 ) I, where θ is the current temperature, θ0 is the initial temperature, α is the linear thermal expansion of material, and I is the identity tensor. The thermal distribution can be achieved by the wellknown energy balance equation. Considering the energy of volumetric strain in the transient energy equation, the fully coupled governing equations of thermo-elasticity can be written for an isotropic domain as ρ c p θ˙ + β θ0 tr (˙ε) = −∇ · q + G h

(3)

where c p is the specific heat parameter at constant pressure, β is a thermo-elastic constant defined as β = α(3λ + 2μ) with λ and μ denoting the Lamé parameters, G h is the rate of internal energy as a source or sink, q is the heat flux vector that is related to the temperature by a proper constitutive law for heat transfer process. In above, θ˙ represents time derivative of θ and tr() is the well-known trace operator of an arbitrary tensor field . The gap between the two faces of discontinuity can be considered as a distinct medium, so the interaction between this medium and the bulk domain is accomplished by transaction of the heat flux through the discontinuous faces. As a result, the heat flux can be considered as the internal boundary condition for the thermal phase. In many applications, the medium between these discontinuous faces contains very low conductive material, such as the air, and the heat flux can be ignored through the discontinuous faces. Hence, the partial differential equation can be associated with the external and internal thermal boundary conditions as

123

Computational Mechanics

q .n = −qn θ = θ˜

on

q

on

θ

q .nd = −qcont on

(4)

c ⊂ d

where qn is the prescribed heat flux on q , θ˜ is the prescribed temperature on θ , and qcont is the heat flux imposed on the discontinuous contact surface d . Note that the heat flux can be transferred only on the contact interface c . In the case of discontinuous opening, the adiabatic boundary condition is imposed on discontinuous surfaces, which results in the free flux condition. The constitutive equation relating the flux vector to temperature can be expressed through the heat transfer analysis. The movement of particles is neglected through the infinitesimal condition of thermo-elasticity and the conduction mode is only considered in the heat transfer process. The conductive heat transfer rule is governed by the Fourier law of conduction defined as q = −κ ∇θ , where κ is the second order conductivity tensor. For the isotropic material, the conductivity tensor can be represented by a scalar constant κ.

2.1 The weak form of governing equations The weak form of thermo-mechanical governing equations can be obtained by defining an appropriate space of test and trial functions corresponding to unknown field variables. The space of the trial function U and corresponding test function δU for the displacement field can be defined as U = {u ∈ S u |u = u¯ on u and u discontinuous on d } (5) δU = {δu ∈ Su |δu = 0 on u and δu discontinuous on d }

(6) where the space Su is related to the regularity of the solution that basically allows for discontinuous functions across the discontinuity [42]. Multiplying the equilibrium equation by the test function, integrating over the domain of the problem, and employing the Divergence theorem for the domain with internal boundary as a discontinuity, the weak form of the governing boundary value problem can be obtained as [43] 

 (δε)T σ d + δuT tcont d  c   T = (δu) ρ b d + (δu)T t d 

t

(7)

in which the conventional Voigt notation is used for derivation of the weak formulation. The above integral equation must be satisfied for all arbitrary test functions δu. The notation δu denotes the jump value of δu across the discontinuity d . It is noteworthy to highlight that with the

123

discontinuous approach, the second integral term is appeared to impose the tractions on the faces of c ⊂ d due to the mechanical contact condition. In a similar manner, the space of the trial function T and corresponding test function δT for the temperature field can be defined as T = {θ ∈ Sθ |θ = θ¯ on θ and θ discontinuous on d

 (8)

δT = {δθ ∈ Sθ |δθ = 0 on θ and δθ discontinuous on d } (9) in which the space Sθ is related to the regularity of the solution that basically allows for discontinuous functions across the discontinuity. Finally, multiplying the energy equation by the test function, integrating over the domain of the problem, and utilizing the Divergence theorem, the weak form of the energy equation can be obtained as follows   δθ T ρc p θ˙ d+ δθ T qcont d+ δθ T βθ0 mT ε˙ d  c     T T = ∇(δθ ) . q d + δθ qn d + δθ T G h d





q



(10) in which the above integral equation must also be satisfied for all arbitrary test functions δθ . In Eq. (10), m is a unit vector defined as mT = [1 , 1 , 0] for 2D and mT = [1 , 1 , 1 , 0 , 0 , 0] for 3D problems. It must be emphasized that the second integral term is appeared to impose the heat flux across the contact interface c ⊂ d due to thermal contact condition at the contact zone. Note that in the derivation of Eq. (10), the free flux condition or the adiabatic boundary condition is imposed on discontinuous opening surfaces.

3 Modeling the thermo-mechanical contact problem The thermo-mechanical contact problem must be characterized by two contact constraints; the mechanical contact constraint and the thermal contact constraint. The former prevents the penetration between two contact surfaces; once the contact condition is recognized the heat transfer occurs across the contact interface. The latter describes the interchangeable heat flux that can be expressed by employing a pressure dependent model, and leads to an implicit coupled relation between the thermal and mechanical contact conditions.

Computational Mechanics

3.1 Frictional contact along the discontinuous surface In order to represent a frictional contact model in the framework of an X-FEM formulation along the discontinuous interface, a vector function g(x) = u(x) is introduced corresponding to the relative displacement between two contact interfaces, as shown in Fig. 2. The normal and tangential gap functions can be defined as g N (x) = χ N (x) . g(x) = g N (x) nd

x ∈ d

gT (x) = g(x) − g N (x) = gT (x) md

(11)

where χ N is the normal projection tensor defined as χ N = nd ⊗ nd and, nd and md denote normal and tangential unit vectors to the discontinuous surface, respectively. The contact constraints can be imposed in the normal direction at the discontinuous surface d based on the Karush-Kuhn-Tucker condition to prevent the penetration of two bodies. If the contact condition occurs the compressive contact forces are induced to satisfy the zero normal gap, i.e. t Ncont = tcont . nd < 0 and g N = 0. Otherwise, the traction free condition is fulfilled due to opening condition, i.e. t Ncont = 0 and g N > 0. These inequality conditions can be summarized using the well-known Karush-Kuhn-Tucker relation as g N ≥ 0 , t Ncont ≤ 0 , g N t Ncont = 0

(12)

In the tangential direction, two conditions of the ‘stick’ and ‘slip’ are distinguished based on the Coulomb’s friction law, defined in two-dimensional space as     F f = tTcont  − μ f t Ncont 



= 0 Slip < 0 Stick

(13)

where μ f is the friction coefficient and tTcont is the tangential contact traction defined as tTcont = tTcont md , with  cont  cont md = tT /tT  denoting the unit tangential vector. For the stick condition the tangential gap function vanishes, i.e. gT = 0 and F f < 0. In the slip condition, the incremental tangential gap is positive, i.e. dgT > 0 and F f = 0.

3.2 Thermal contact along the discontinuous surface There are basically two types of contact conditions employed in the context of thermal conduction; a perfect and an imperfect contact interface. If the temperature and normal heat flux are continuous across the interface, it is called as the perfect interface, or the iso-thermal contact condition; otherwise, the interface is imperfect. Since the assumptions of perfect interface may not be valid in many practical problems, an imperfect contact interface is required. According to the theory of thermal imperfect interfaces, there are different thermal contact conditions, which are known as the Kapitza’s thermal interface, the highly conducting imperfect interface, and the general imperfect interface [44]. The widely adopted Kapitza’s thermal interface model affirmed that the temperature has a jump across the interface but the normal heat flux is continuous and proportional to the temperature jump [45– 47]. In highly conducting interface model, the temperature is continuous while the normal heat flux is discontinuous across the interface; in such case, the jump of heat flux discontinuity is proportional to the temperature gradient of the surface and can be derived using the Laplace–Young equation [22,48]. In the general imperfect interface, the temperature and normal heat flux are both discontinuous across the interface. Consider the low conductive condition, e.g. air, between two discontinuous surfaces, the heat transfer in the thermal contact zone can be described by the Kapitza’s constitutive law as qcont (x) =

1

θ (x) Rcont = h cont θ (x)

x ∈ c

(14)

(Slave)

x

d

gN ( x)

nd

g( x ) g T ( x)

x

d

(Master) Fig. 2 Description of the frictional contact within the X-FEM framework along the discontinuous interface

where Rcont and h cont are the thermal resistance and thermal conductivity coefficients of the interface, respectively. In a typical extreme condition where the thermal conductivity is close to zero, the well-known adiabatic assumption can be assumed. The conductivity coefficient can be affected by different parameters, such as the contact pressure, material property, temperature, etc. A micro-mechanical study is generally used to determine the effect of conductivity coefficient in the macroscopic level [49]. At the microscopic level, the real contact surface is assumed as a rough surface, where the contact takes place at certain spots, and as a result the real contact zone of microscopic level is less than the apparent contact area in macroscopic level. Hence, the heat exchange

123

Computational Mechanics

between two contacting surfaces can be occurred in the following three types; the direct conduction through the spots, the conduction via the air contained in the cavities, and the radiation between cavity surfaces. Since the air in cavities has a low conductive condition and the radiation between the cavity surfaces needs a very high variation of the temperature, the spots can be taken as the most significant manner for the exchange of heat between two surfaces. If the contact zone consists of Nsp spots with the conductivity of h i for each spot, the total conductivity of spots can be obtained using the assumption of parallel effect of spots at the con Nsp tact zone by h cont = i=1 h i [49]. Applying an equivalent spot in the microscopic level, the total conductivity can be ˆ with hˆ denoting the conductivity of defined by h cont = Nsp h, an equivalent spot. There are various experimental and analytical relations reported in the literature for determination ˆ From the analytical point of view, the solution of the of h. conduction equation can be used in a specific domain of the ˆ One of the most practical spot area to derive a relation for h. relations employed for description of the heat flux across a contact area is defined as [32]

the aim of this study is focused on the quasi-static formulation of continuum mechanics in an infinitesimal strain framework, as shown through numerical examples in Sect. 7. Thus, the sliding velocity is neglected and the generated heat flux qfrict is removed from the total energy Q total cont = qcont + qfrict over the contact surface.

 h cont t Ncont (x, t)

cont ε t (x, t) c = h c0 N He

in which λ N > 0 is also known as the normal penalty parameter. For the tangential contact constraint, the tangential contact stiffness λT is imposed in an incremental manner using the return mapping algorithm within the framework of frictional theory of plasticity. In this case, the stick condition is enforced by applying a large value of the tangential penalty parameter λT to prevent the tangential sliding, and the slip condition is imposed by the corresponding value of tangential penalty parameter λT calculated from the return mapping algorithm. Hence, the incremental tangential traction can be computed using the tangential contact stiffness λT as

x ∈ c

(15)

where h c0 , He and εc are respectively the resistance coefficient, the Vickers hardness, and the material constant. Obviously, a higher contact pressure, i.e. a larger value of t Ncont , and a soften material, i.e. a smaller value of He , provide a larger contact area and as a result, a greater heat conductivity. Note that the thermal contact condition is defined at the opening surfaces where the contact heat flux is vanished, i.e. g N > 0 and qcont = 0, while in the case of contact condition, i.e. g N = 0, the thermal contact condition is defined by relations (14) and (15) to describe the heat flux across the contact surface. Thus, the thermal contact condition can be defined as g N ≥ 0 and qcont .g N = 0. It is noteworthy to highlight that there is also another source of thermal energy at the contact surface due to frictional sliding. In fact, frictional sliding transforms the dissipated elastic energy to thermal energy over the contact surface. The generated heat flux qfrict can be defined as   qfrict (x) = tTcont (x) g˙ T (x)   = μ f t cont (x) g˙ T (x) N

x ∈ c

(16)

where g˙ T is the slipping velocity [35]. In some applications, it is reasonable to ignore the effect of thermal energy due to frictional sliding qfrict in comparison with the heat transmitted between two contact faces qcont , where the slipping velocity is not significant, such as rock problems [50–52]. Moreover,

123

3.3 Imposition of contact constraints In order to impose the Karush–Kuhn–Tucker condition, i.e. g N ≥ 0, t Ncont ≤ 0, g N .t Ncont = 0, together with the stick/slip contact condition, i.e. dgT ≥ 0, F f ≤ 0, dgT .F f = 0, and the thermal contact condition, i.e. g N ≥ 0, qcont .g N = 0, the penalty method is employed here for each contact constraint. The normal contact constraint is enforced to prevent the penetration between two contact surfaces. Hence, the normal contact traction can be computed employing the normal contact stiffness λ N as t Ncont (x, t) = λ N .g N (x, t);

lim g N (x, t) = 0

λN → ∞

dtTcont (x, t) = λT (x, t). dgT (x, t);  λT (x, t)|F f ≤ 0 , λT > 0

(17)

(18)

in which the elastic-predictor/plastic-corrector is employed using a proper value of the tangential penalty parameter λT to impose the stick/slip contact condition. For each increment, the elastic (stick) predictor is first assumed employing a large value of λT , the plastic corrector is then applied to distinguish between the stick and slip conditions. If the slip condition occurs the tangential penalty parameter λT is modified as   μ f t Ncont (x, t) − tTcont (x, t) (19) λT (x, t) = dgT (x, t) Note that the proper values of the normal and tangential (stick) penalty parameters play a significant role to avoid the singular solution of the governing equations. Finally, the thermal contact constraint, i.e. g N ≥ 0, qcont .g N = 0, is employed to physically implement a virtual thermalconductor with the nonlinear thermal conductivity (15). In

Computational Mechanics

fact, relation (14) represents a penalty function where the nonlinear penalty parameter is defined in (15) as the nonlinear thermal conductivity. Hence, the thermal contact constraint can be summarized as  qcont (x, t) = h cont t Ncont (x, t) .θ (x, t) with

cont ε t N (x, t) c h cont (x, t) = h c0 He

(20)

4 X-FEM formulation of thermo-mechanical contact problem In order to numerically model the thermo-mechanical contact problem described in preceding sections, the extended finite element method is employed in which the contact interface is embedded within the finite element mesh by adding appropriate enrichment functions for the mechanical and thermal approximations. The extended finite element method is employed to discretize the weak form of governing Eqs. (7) and (10) in spatial domain. In X-FEM, the basis functions of mechanical and thermal approximations are enriched by special discontinuous functions through the partition of unity method, where additional degrees of freedom are used to model the discontinuity embedded within the finite element mesh. The spatially discretized governing equations are then integrated in time using the Newmark time splitting method. The final systems of nonlinear algebraic equations are then solved using the Newton–Raphson iterative scheme in a staggered and fully coupled manners.

4.1 Approximation of displacement field In order to incorporate the displacement jump across the contact interface, the trial and test functions of the displacement field, i.e. u(x, t) and δu(x, t), are respectively approximated using the extended finite element discretization. The Heaviside enrichment function is used to model the discontinuous displacement field. Hence, the X-FEM approximation of the displacement field uh (x, t) is defined based on an additive decomposition of the approximation into the standard and enriched parts as uh (x, t) =

Nui (x)u¯ i (t)

+

H (ϕ(x)) =

+1

ϕ(x) ≥ 0

0

ϕ(x) < 0

(22)

where ϕ(x) is the signed distance function, which is represented by the level set function in the context of X-FEM. The signed distance function is defined based on the absolute distance value of a point x from the contact interface as ϕ(x) = x − xd sign ((x − xd ) · nd (xd ))

(23)

where xd is a point on the contact interface d which has the closest distance from the point x, and nd is the normal vector to the contact surface at point xd . Note that the displacement jump across the discontinuity d can be expressed by uh (x, t) = Nu (x).¯a(t). The strain vector εh (x, t) corresponding to the approximate displacement field can be written in terms of the standard and enriched nodal values by differentiating from the enriched displacement field (21) as εh (x, t) = Bu (x) u(t) ¯ + Ba (x)¯a(t)

(24)

where Bu (x) = SNu (x) and Ba (x) = SNa (x) are respectively the spatial derivatives of the standard and enriched shape functions, with S denoting the strain differential operator.

4.2 Approximation of temperature field In order to arrive at the finite element approximation of the temperature field θ h (x, t), the Heaviside enrichment function is also employed to model the temperature jump across the contact interface within the X-FEM framework. Hence, the X-FEM approximation of the temperature field θ h (x, t) can be defined as

Nθi (x)θ¯i (t)

i∈ N std

Nui (x) (H (ϕ(x)) − H (ϕ(xi ))) a¯ i (t)

+

i∈N enr

¯ + Na (x)¯a(t) = Nu (x) u(t)



θ h (x, t) =

i∈N std

are the standard displacement DOFs, and a¯ i (t) are additional displacement DOFs corresponding to the Heaviside enrichment function. Because of independent displacement fields at the two sides of the contact surface, the Heaviside jump function is an appropriate enrichment function for simulation of the contact problem defined as

(21)

where Nui (x) is the standard displacement shape function of node i, N std is the set of standard nodal points, N enr is the set

of nodes whose support is cut by the contact interface, u¯ i (t)

i∈ N enr

Nθi (x) (H(ϕ(x)) − H(ϕ(xi ))) ω¯ i (t)

¯ = Nθ (x)θ¯ (t) + Nω (x)ω(t)

(25)

where Nθi (x) is the standard temperature shape function of node i, θ¯i (t) are the standard nodal temperature DOFs, and

123

Computational Mechanics

ω¯ i (t) are additional nodal temperature DOFs corresponding to the Heaviside enrichment function. It must be noted that the temperature jump across the contact interface can be ¯ expressed by θ h (x, t) = Nθ (x) ω(t). The temperature gradient vector ∇θ h (x, t) can be approximated by taking derivation from the temperature field (25) as ¯ + Bω (x)ω(t) ¯ ∇θ h (x, t) = Bθ (x)θ(t)

4.3 Discretization of governing equations The discrete form of the thermo-mechanical contact problem can be obtained by substituting the approximate displacement and strain fields uh (x, t) and εh (x, t) from (21) and (24) together with the approximate temperature and temperature gradient fields θ h (x, t) and ∇θ h (x, t) from (25) and (26), and their approximate test functions, i.e. δuh (x, t), δεh (x, t), δθ h (x, t) and δ∇θ h (x, t), into the weak form of governing Eqs. (7) and (10), and satisfying the necessity that the weak form should hold for all arbitrary test functions. According to the Bubnov–Galerkin technique, the approximate test functions δuh (x, t) and δε h (x, t) are considered in the same approximating space of uh (x, t) and ε h (x, t) as ¯ + Na (x)δ a¯ (t) δuh (x, t) = Nu (x) δ u(t)

(27a)

¯ + Ba (x)δ a¯ (t) δε (x, t) = Bu (x) δ u(t)

(27b)

h

δθ h (x, t)

Similarly the approximate test functions and δ∇θ h (x, t) are considered in the same approximating space of θ h (x, t) and ∇θ h (x, t) as ¯ δθ h (x, t) = Nθ (x)δ θ¯ (t) + Nω (x)δ ω(t) h ¯ ¯ δ∇θ (x, t) = Bθ (x)δ θ(t) + Bω (x)δ ω(t)

(28a) (28b)

Substituting the enriched FE approximations of the test and trial functions of displacement and strain fields into the weak form of governing Eq. (7), the discretized mechanical governing equations can be obtained as =

1 2



 =

ext thm Fint u − Fu − Fu ext thm cont Fint a − Fa − Fa + F

 Faext

=

 =0

(29)

Fuint

=

Faint

=

(Bu ) σ d T

 

(Ba )T σ d

123



t

(Na ) ρ b d +

t

 Futhm

=



Fathm = F

(Na )T t d

(Bu )T σθ d



(Ba )T σθ d



=

cont

(Nu )T t d





c

(Nu )T tcont d

(30)

in which the internal force vector F cont is defined over the contact interface c that takes into account the effect of internal forces emanated from the contact condition between two discontinuity faces. Moreover, the force vectors Fthm represent the contribution of thermal stress σθ on the mechanical deformation, i.e. σθ = Dεθ . In a similar manner, substitute the enriched FE approximations of the test and trial functions of the temperature and temperature gradient fields into the weak form of governing Eq. (10), the discretized thermal governing equations can be obtained as   1 = 2   ext ˙ strg ˙ mech Qint θ + Qθ − Qθ − Qθ = = 0 (31) ext ˙ strg ˙ mech − Qcont Qint ω + Qω − Qω − Qω where  Qint θ = Qint ω



=



strg Qθ

(Bω )T q d



Qext θ = Qext ω

(Bθ )T q d

 

=



(Nθ )T G h d +



q

(Nω ) G h d +



 =

strg



Qω = Qmech θ

=

Qmech ω

=

Qcont =



 

q

(Nω )T qn d

(Nθ )T qstrg d (Nω )T qstrg d



 

(Nθ )T qn d

 T

where 

 (Nu )T ρ b d + T

(26)

where Bθ (x) = LNθ (x) and Bω (x) = LNω (x) are respectively the gradient of the standard and enriched temperature shape functions, and L is the well known gradient operator.



 Fuext =



c

(Nθ )T qmech d (Nω )T qmech d (Nθ )T qcont d

(32)

in which the internal flux vector Qcont is defined over the surface integral c that takes into account the effect of inter-

Computational Mechanics

nal fluxes transferred through the contact surface between two discontinuity faces. Moreover, the flux vectors Qmech and Qstrg are defined respectively to represent the contribution of deformation rate in the mechanical energy, i.e. qmech = βθ0 tr (˙ε), and the heat capacity in the storage energy, i.e. qstrg = ρ c p θ . In order to discretize the X-FEM Eqs. (29) and (31) in time, the Newmark time splitting method is employed using the GN11 time integration scheme for the velocity field u˙ and the time derivative of temperature field θ˙ . The link between the successive values of the unknown field variables at time tn+1 = tn + t and the known field variables at time tn can be established as 

1 1 u˙ n+1 = − 1 u˙ n − un ) − (33a) (u γ t n+1 γ 

1 1 θ˙n+1 = − 1 θ˙n − θn ) − (33b) (θ γ t n+1 γ where 0 < γ < 1 is the Newmark parameter, in which for unconditionally stability of the time splitting algorithm, it is required that γ ≥ 0.5. Substituting relations (33) into Eqs. (29) and (31) results in   1  n+1 =  2 n+1   Fuint − Fuext − Futhm = =0 (34) Faint − Faext − Fathm + F cont n+1 

⎡ ⎢  J =⎢ ⎣ 

BuT D Bu

d

BaT D Bu d

 

n+1 =

1 2



 =

k k k  k+1 n+1   n+1 + Jn+1 dXn = 0

 ∂ 1 J =

d



BaT D Ba d + Kcont



ext Qint θ + Qθ

ω

ω

(38)

d

 

BuT βm Bω

⎤

d ⎥ ⎥ ⎦ T Ba βm Bω d

(39)

where the contact stiffness matrix Kcont is defined as



Kcont = =

∂F cont  ∂ a¯ c



n

(35) ¯ is defined as where Q     strg ¯θ Q Qθ + Qmech 1 θ = γ t Qstrg + Qmech ¯ω Q ω ω   ˙ strg ˙ mech 

Qθ + Qθ 1 −1 + γ ˙ mech ˙ strg + Q Q



or

BaT βm Bθ d

cont ext Qint ω + Qω − Q n+1 n+1     strg mech ¯θ Q Qθ + Qθ 1 − + =0 strg γ t Q + Qmech ¯ω Q ω ω n+1

∂ 1 ∂ 1 ∂ 1 ¯ ∂ u¯ ∂ a¯ ∂ θ¯ ∂ ω ∂ 2 ∂ 2 ∂ 2 ∂ 2 ¯ ∂ u¯ ∂ a¯ ∂ θ¯ ∂ ω

 BuT βm Bθ

(37)

where J is the Jacobian matrix, and dX is the vector of incremental unknown nodal values defined as dXT =   T T T T ¯ . The solution of the linearized system d u¯ d a¯ d θ¯ d ω of equations (37) results in an increment of the standard and enriched nodal DOFs, in which the corresponding unknown k+1 k = Xn+1 + dXnk . The nodal values can be obtained by Xn+1 k k Jacobian matrix can be derived from ∂ n+1 /∂Xn+1 as

 BuT D Ba

and 

The system of discretized Eqs. (34) and (35) are nonlinear because of the frictional contact behavior and nonlinear thermal contact condition. However, it must be emphasized that even for the frictionless contact condition and linear thermal conductivity of contact surface, the problem is still kinematically nonlinear, since the active contact zone c is not known a priori. As a consequence, the Newton–Raphson iterative procedure is employed to linearize the system of Eqs. (34) and (35) at each iteration k of time interval (tn , tn+1 ). By expanding the residual Eq. (34) with the first-order truncated Taylor series, the linearized form of nonlinear system of Eq. (34) can be obtained as

(36)

=

c

(Nu )T

∂  cont t N nd + tTcont md Nu d ∂ a¯

(Nu )T Dcont Nu d

(40)

in which the contact property matrix is defined as Dcont = λ N nd ⊗ nd + λT md ⊗ md . Similarly, the residual Eq. (35) can be expanded by the first-order truncated Taylor series to obtain the linearized system of equations as k k ¯k k+1 n+1  n+1 + J n+1 dXn = 0

(41)

123

Computational Mechanics

where the Jacobian matrix J¯ can be derived from ∂kn+1 / k as ∂Xn+1 ⎡ J¯ = ⎣

∂1 ∂1 ∂1 ∂1 ¯ ∂ u¯ ∂ a¯ ∂ω ∂ θ¯ ∂2 ∂2 ∂2 ∂2 ¯ ∂ u¯ ∂ a¯ ∂ω ∂ θ¯

⎤

Thus, the coupling matrix Hcoup in (45) can be defined as H

⎦



 coup

(42)

=

c

(Nθ )

T



θ 

h c0 εc He



t Ncont He

εc −1 

(nd )T Nu d

 λN (48)

or 

⎡ ⎢ ⎢ J¯ = ⎢ ⎣

1 γ t 1 γ t







 NθT βm Bu d

1 γ t

NωT βm Bu

1 γ t

1 γ t





NθT βm Ba d NωT βm Ba d + Hcoup



 1 γ t

d



NθT ρcNθ



d +

NωT ρcNθ d +



 



 BθT κBθ

d

1 γ t

BωT κBθ d

1 γ t

 

NθT ρcNω

d +

NωT ρcNω d +

where the thermal contact conductive matrix Hcont is defined as ∂Qcont Hcont = ¯  ∂ω ∂qcont d (Nθ )T = ¯ ∂ω   c = (Nθ )T h cont Nθ d (44) c

in which the thermal contact conductivity h cont is defined in relation (15) as a function of the normal contact pressure. In relation (43), the coupling matrix Hcoup represents the effect of thermal contact on the displacement field, which is defined as ∂Qcont Hcoup =  ∂ a¯ ∂qcont d (45) (Nθ )T = ∂ a¯ c By taking the partial derivation from the contact heat flux (20), it results in ∂qcont ∂qcont ∂h cont ∂t Ncont = ∂ a¯ ∂h cont ∂t Ncont ∂ a¯ 

ε −1  cont ∂t N h c0 εc t Ncont c = θ  He He ∂ a¯

(46)

Now, taking the derivation from relation (17), it leads to a relationship between the contact pressure and the displacement normal gap as ∂t Ncont ∂ gN = λN = λ N (nd )T Nu ∂ a¯ ∂ a¯

123

(47)

 

BθT κBω

⎤ d

BωT κBω d + Hcont

⎥ ⎥ ⎦

(43)

5 The staggered Newton solution algorithm The set of Eqs. (37) and (41) must be solved simultaneously to obtain the vector of incremental unknown nodal values dX. In this manner, the coupled system of nonlinear equations can be written as ⎫k ⎧ ⎡ J uu 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎬ ⎨ 2 ⎪ ⎢ J au +⎢ ⎢ J¯ ⎪ ⎪ 1 ⎪ ⎪ ⎣ θu ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ 2 n+1 J¯ ωu

⎫k ⎧ d u¯ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎥ ⎪ ⎬ ⎨ d a¯ ⎪ J aω ⎥ ⎥ =0 ⎪ d θ¯ ⎪ J¯ θω ⎥ ⎪ ⎪ ⎦ ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ ¯ n J¯ ωω n+1 d ω

J ua J uθ J uω J aa J aθ J¯ θa J¯ θθ J¯ ωa J¯ ωθ

⎤k

(49) In order to solve the nonlinear coupled Eq. (49), two techniques can be employed; a monolithic approach, in which the coupled equations are solved simultaneously, and a staggered approach, in which a partitioned solution strategy is used to solve the mechanical and thermal equations individually. Basically, the staggered approach has several advantages over the monolithic approach from computational point of view; firstly, each set of governing equations can be solved by particular solvers, so the need for a unified direct solver that is computationally expensive can be eliminated. Secondly, the staggered algorithm offers desirable flexibility from a computer programming viewpoint; as a result, the thermal governing equations can be easily implemented within an existing mechanical FE computer code. Thus, a staggered Newton solution strategy is employed here to solve the set of nonlinear coupled Eqs. (37) and (41) in an iterative sequential algorithm. In the staggered Newton solution algorithm, the nonlinear thermal Eq. (41) are first solved at each iteration k of time interval (tn , tn+1 ) to obtain the incremental temperature field

Computational Mechanics

by assuming an incremental displacement field; the nonlinear mechanical Eq. (37) is then solved by updating the nodal temperature values to obtain the nodal displacement values. This procedure is continued at each time step until the convergence is achieved within a predetermined error tolerance. Note that the stability and convergence of the solution can be achieved through the proposed staggered Newton algorithm [53–55]. It must be emphasized that the linearization of nonlinear mechanical and thermal equilibrium equations allows to directly impose the inequality contact constraints within each iteration, i.e. the Karush-Kuhn-Tucker condition g N ≥ 0, t Ncont ≤ 0, and g N .t Ncont = 0 to prevent the penetration between two contact surfaces, the stick/slip contact condition dgT ≥ 0, F f ≤ 0, and dgT .F f = 0 to capture the stick/slip behavior between two contact surfaces in tangential direction, and the thermal contact condition g N ≥ 0, and qcont .g N = 0 to impose the nonlinear thermal conduction between two contact surfaces, as demonstrated in Sect. 3.3. According to the procedure described above, the nonlinear coupled Eq. (49) are first solved for the thermal phase to obtain the incremental temperature field by assuming an incremental displacement field as k  k  d θ¯ J¯ θθ J¯ θω ¯ n J¯ ωθ J¯ ωω n+1 d ω   k  k  k J¯ θu J¯ θa 1 d u¯ =− − (50a) ¯ ¯ 2 d a¯ J ωu J ωa 

θ¯ ¯ ω

n+1



k+1

= n+1

θ¯ ¯ ω

k



+ n+1

d θ¯ ¯ dω

k

n+1

n

(50b) n

In the next step, the nonlinear coupled Eq. (49) are solved for mechanical deformations by updating the nodal temperature values to obtain the updated nodal displacement values as k  k  d u¯ J uu J ua J au J aa n+1 d a¯ n   k k  k 1 J uθ J uω d θ¯ =− − ¯ n  2 n+1 J aθ J aω n+1 d ω  k+1  k  k u¯ u¯ d u¯ = + a¯ n+1 a¯ n+1 d a¯ n

(51a)

(51b)

Finally, the iterative sequential procedure is continued until the convergence is achieved within an acceptable predetermined error. Note that in the numerical simulation examples presented in Sect. 7, the effect of mechanical deformation in producing the thermal energy is neglected because of the low thermal expansion coefficients and low loading rate proposed here [18–21,56].

6 Numerical integration algorithm One of the main issues with the X-FEM technique is the numerical integration of an element in the presence of discontinuity in the displacement and temperature fields. The discontinuous nature of the solution in this case requires a special treatment to accurately integrate the non-smooth enrichment functions appearing in the enriched part of the approximation. There are basically two approaches for numerical integration in the X-FEM framework; the rectangular sub-girds method and the triangular partitioning method (Fig. 3). Implementation of the rectangular sub-girds method is more straightforward than the triangular partitioning method, however—the accuracy of latter is more reliable than the former. This is due to the fact that in the triangular partitioning method the numerical integration of an element bisected by the discontinuity is performed separately on each sub-element into which the element is divided. In this approach, the integration over each sub-element on either side of the discontinuity is partitioned into triangles over which the standard Gauss quadrature is applied. In this manner, the numerical integration is employed by an isoparametric mapping of the Gauss quadrature points into the sub-elements, as shown in Fig. 3. The remaining elements which are not cut by the discontinuity are integrated using the conventional Gauss quadrature scheme. In the next section, the numerical integration of the first four examples is performed using four rectangular sub-girds for elements bisected by the discontinuity, in which four Gauss quadrature points are used in each sub-grid. However in the last example, the triangular partitioning method is employed by partitioning each sub-element into five triangles, in which three Gauss quadrature points are used in each sub-triangle. Moreover, the numerical integrations of mechanical contact Eq. (40) and thermal contact Eqs. (44) and (45) are significant and must be computed properly. To this end, a number of Gauss points are set along the one-dimensional segments of the discontinuity delimited with the element edges to integrate the mechanical and thermal contact terms along the discontinuity. For the integration of mechanical contact Eq. (40), it is obvious that the two Gauss quadrature points can be used adequately to evaluate the numerical integration of the discontinuity since the bilinear basis functions are employed for the standard degrees of freedom [10,11,16]. So, the numerical integration is performed with two quadrature points over each segment of the fracture discontinuity which is cut by an element, as shown in Fig. 3a. However for numerical integration of the thermal contact Eqs. (44) and (45), it is impossible to indicate the required number of quadrature points, since the order of conductivity function (15) is not determined in prior and is dependent on εc . For this reason, a sufficient fine mesh is employed together with two Gauss quadrature points to compute the numerical integra-

123

Computational Mechanics

η

Jξ

ce

(a)

ξ

e

Line integration points

e

ξ′

Integration points of sub-grid

Partitioning of mapped domain

Jξ

Γce

η′

Jξ ′

η

Jξ ′

η′

ξ

ξ′ (b) Fig. 3 Numerical integration of an enriched element along the fracture interface and over the element using a the rectangular sub-girds method and b the triangular partitioning method

tion of thermal contact equations, since the thermal contact constitute law is a smooth and piecewise linear function. Thus, it is important to appropriately evaluate the thermal contact integrations by choosing a sufficiently fine mesh around the discontinuity, or increasing the number of Gauss quadrature points along the fracture discontinuity.

7 Numerical simulation results In order to illustrate the accuracy and robustness of the proposed computational algorithm, several thermo-mechanical contact problems are solved numerically, including: the numerical simulation of thermal resistance in a fractured body with constant conductivity at the contact surface, thermo-mechanical modeling of a fractured body with a pressure-dependent thermal contact model, thermomechanical modeling of a fractured body with frictional contact behavior, thermo-mechanical modeling of a doubleclamped beam with a central crack, and numerical modeling of a complex fractured geometry for a plate with multiple cracks subjected to thermo-mechanical loading conditions. All these examples are chosen to verify the accuracy of the proposed X-FEM thermo-mechanical contact algorithm for modeling the thermal resistance in a fractured body, the mechanical behavior in a thermal contact problem, and the combined thermo-mechanical contact behavior in friction and frictionless contact problems. The numerical results are

123

compared with those reported in the literature, or available analytical solutions. It is noteworthy to mention that at each Newton–Raphson iteration, the contact stiffness matrix Kcont , the contact force vector F cont , the thermal contact conductive matrix Hcont , and the contact flux vector Qcont associated with each segment are incorporated into the final system of equations if the set of integration points over the contact segment has a negative normal jump, i.e. g N < 0. Moreover, although the unconditionally stable condition of the Newmark method is generally proposed by choosing γ ≥ 0.5, the accuracy of the method significantly depends on the proper choice of the Newmark parameter and the time step size. In this study, the Newmark parameter is chosen as γ = 0.65, which is very close to that previously suggested by Zienkiewicz et al. [57] in the dynamic FEM analysis and Khoei et al. [19] within the dynamic X-FEM framework. In the examples proposed in this section, it is shown that this value leads to a better accuracy of the solution. On the other hand, the proper accuracy of the solution can be obtained by choosing the Newmark parameter close to the lower bound of unconditionally stability and selecting an appropriate time increment.

7.1 Modeling the thermal resistance in a fractured body The first example is chosen to verify the accuracy of the proposed computational algorithm in simulation of the thermal

Computational Mechanics 1

fracture interface 0.8

T = 1o C

T =0 κ = 1.7 W mC c = 880 J kgC

y

ρ = 2650 kg m3

x Fig. 4 Modeling the thermal resistance in a fractured body with constant conductivity at the contact surface; the geometry, boundary conditions, and material properties

Temperature (°C)

hcont

0.6

hcont = 1.0 E + 5 hcont = 5.0 hcont = 0.2

0.4

hcont = 1.0 E − 5

0.2

0 0

contact phase by modeling the thermal resistance in a fractured body with constant thermal conductivity at the contact zone. This problem was originally proposed by Yvonnet et al. [23] for modeling the Kapitza thermal resistance within the X-FEM. The geometry, boundary conditions, and material properties are presented in Fig. 4. The temperature is set to zero at the left edge and T = 1 ◦ C at the right edge of the domain. The initial temperature is set to zero in the whole domain. A uniform structured FE mesh is utilized with 19 × 19 quadrilateral elements. A vertical discontinuity is assumed at the middle of the body, and the fully contact condition is employed with a constant thermal conductivity h cont at the contact interface. In Fig. 5, the evolutions of temperature are plotted for various thermal conductivities at two sampling points located in the center of left and right blocks. Obviously, at the beginning of simulation, the contact interface does not affect the temperature evolutions of two blocks for various thermal conductivities, however the effect of thermal contact conductivity becomes significant when the time increases and the heat transfer occurs across the contact interface. It can be observed that by increasing the thermal contact conductivity the temperature difference between two sampling points becomes closer. Clearly, it is obvious that the higher value of thermal conductivity leads to increase the heat flux across the contact interface from the right block to left one. Obviously, the thermal contact conductivity has a significant impact on the temperature difference between two sampling points in the range of 0.2 < h cont < 5. In Fig. 6, the evolutions of temperature are plotted along the horizontal direction for various thermal conductivities, i.e. h cont = 1.0E + 5, 5.0, 0.2, 1.0E − 5. A comparison is performed in this figure between different transient dynamic analyses and the steady state solution. It can be seen from the plots of temperature field that there is a jump across the contact interface as the thermal contact conductivity of fracture decreases, particularly a fully continuous temperature

0.5

1

1.5

Time

2

2.5

3

(×106s)

Fig. 5 The evolutions of temperature for various thermal conductivities h cont at two sampling points in the center of left and right blocks; the dashed line belongs to the point at the center of left block, and the solid line belongs to the point at the center of right block

field with h cont = 1.0E + 5, presented in Fig. 6a, transfers to a completely discontinuous temperature field with h cont = 1.0E − 5, displayed in Fig. 6d. Obviously, the Heaviside enrichment function is used efficiently in X-FEM modeling of the temperature field across the contact surface for both the continuous temperature field with a large value of the thermal contact conductivity and the discontinuous temperature field in the case of a low value of thermal contact conductivity. Furthermore, the results of transient dynamic analyses are compared in Fig. 6 with those of analytical steady state solutions. The analytical solution of temperature field can be straightforwardly derived as T (x) =

⎧ ⎨ ⎩

κ R (TR −TL )h cont (κ L +κ R )L h cont+κ L κ R (x)+TL κ L (TR −TL )h cont (κ L+κ R )L h cont+κ L κ R (x − 2L)+TR

for 0 < x < L for L < x < 2L

(52) where TL and TR are respectively the left and right boundary temperatures of the body, L is the length of each block, and κ L and κ R are correspondingly the thermal conductivities of the left and right blocks. Obviously, a complete agreement can be seen between the proposed numerical simulations and those of analytical steady state solutions.

7.2 Thermo-mechanical analysis of a fractured body with the pressure-dependent thermal contact model The second example is chosen to illustrate the robustness of the proposed computational algorithm in the thermo-

123

Computational Mechanics

(b)

1

t = 0.02 E + 6 s t = 0.2 E + 6 s t = 0.3E + 6 s t = 2.8 E + 6 s

Temperature ( o C)

0.8

1 t = 0.02 E + 6 s t = 0.2 E + 6 s t = 0.3E + 6 s t = 2.8 E + 6 s

0.8

Temperature ( o C)

(a)

steady state

0.6

0.4

steady state

0.6

0.4

0.2

0.2

0

0 0

0.5

1

1.5

0

2

0.5

x − coordinate (m)

(c)

(d)

1 t = 0.02 E + 6 s t = 0.2 E + 6 s

0.6

1.5

2

1.5

2

1 t = 0.02 E + 6 s t = 0.2 E + 6 s t = 0.3E + 6 s

0.8

t = 0.3E + 6 s t = 2.8 E + 6 s steady state

Temperature ( o C)

Temperature ( o C)

0.8

1

x − coordinate (m)

0.4

t = 2.8 E + 6 s steady state

0.6

0.4

0.2

0.2

0

0 0

0.5

1

1.5

2

x − coordinate (m)

0

0.5

1

x − coordinate (m)

Fig. 6 The evolutions of temperature along the horizontal direction for various thermal conductivities; a comparison between different transient dynamic analyses and the analytical steady state (exact) solution; a h cont = 1.0E + 5, b h cont = 5.0, c h cont = 0.2, d h cont = 1.0E − 5

mechanical analysis of a fractured body using a pressuredependent thermal contact conductivity model. This problem was originally proposed by Wriggers and Miehe [32] to present the performance of contact constraints within the coupled thermo-mechanical analysis. The geometry, boundary conditions, and material properties are shown in Fig. 7. In the steady state condition it is sufficient to choose a very large single time step, i.e. t → ∞, and set the Newmark parameter to unity, i.e. γ = 1. The prescribed load p is increased monotonically at each increment of the mechanical phase, and the thermal phase is then solved according to the staggered Newton solution strategy, described in Sect. 5. In order to compare the results of numerical simulation with those of analytical solution, the assumption of stick condition is employed along the contact interface, and the penalty numbers are assumed as η N = ηT = 5.0E +7MN/m3 to enforce the normal and tangential (stick) contact conditions. A uni-

123

form structured FE mesh of 19 × 19 quadrilateral elements is employed for all numerical simulations. In Fig. 8, the variations of temperature field with the contact pressure are plotted at two edges of contact interface together with the evolution of heat flux at the mid-point of the fracture. It is clear that the increase of contact pressure leads to higher value of the thermal contact conductivity, and as a result causes the continuous temperature field across the contact interface. It can be seen from Fig. 8a that a sophisticated pressure-dependent thermal contact conductivity model like the one given in relation (15) is justifiable for a range of the contact pressure less than 1 GPa according to the proposed contact surface properties. On the other hand, the assumption of constant thermal contact conductivity, which leads to extremely reduction of computational costs for such highly nonlinear problem, is only acceptable for the contact pressure value greater than 1 GPa. In order to perform a compari-

Computational Mechanics

p T = 100o C

fracture interface hcont E = 0.07MPa , ν = 0.3

κ = 150 J ms o C

T =0 Fig. 7 The unilateral thermo-mechanical analysis of a fractured body with a pressure-dependent thermal contact model

son between the results of numerical simulation and those obtained from the analytical solution, the problem is solved analytically in the one-dimensional condition. According to Wriggers and Miehe [32], the analytical values of temperature at the top and bottom edges of the contact interface can be obtained as (κ + h cont )Ttop + h cont Tbot κ + 2h cont (κ + h cont )Tbot + h cont Ttop = κ + 2h cont

T+ = T−

(53)

where Ttop and Tbot are respectively the prescribed temperatures at the top and bottom edges of the body and, T + and T − are correspondingly the temperature values at the top

(a)

and bottom edges of the contact interface. Figure 8a presents an excellent agreement between the results of the proposed computational algorithm and those obtained from the analytical solution using a pressure-dependent thermal contact model. It must be noted that coupling between the thermal and mechanical deformations is due to two sources; firstly, the effect of mechanical deformation on the thermal phase due to implementation of a pressure-dependent thermal conductivity on the contact surface, and secondly the influence of thermal effects on the mechanical deformation due to contribution of thermal stresses into the mechanical governing equations, in which the latter can be computed using the thermal expansion coefficient. In order to investigate the effect of thermal expansion coefficient, the evolutions of vertical displacement and normal contact pressure are plotted in Fig. 9 for various thermal expansion coefficients α at different loading steps. Obviously, for low values of the thermal expansion coefficient, i.e. α → 0, the vertical displacement decreases when the external loading increases. However for large values of thermal expansion coefficient, i.e. α = 23.86E − 4 (1/◦ C), the thermal effect becomes dominant in thermo-mechanical deformations of the system for the values of contact pressure between zero and 1 GPa, as shown in Fig. 9a. Moreover, it can be clearly seen from the thermal expansion effect that the tensile thermal stress leads to reduction of the normal contact pressure, as shown in Fig. 9b.

7.3 Thermo-mechanical modeling of a fractured body with frictional contact behavior In the next example, the performance and accuracy of the proposed computational algorithm is presented in thermo-

(b) bottom edge (X-FEM) bottom edge (analytical)

Contact pressure (GPa)

Heat flux ( kN ms )

Temperature ( o C)

top edge (X-FEM) top edge (analytical)

Contact pressure (GPa)

Fig. 8 Thermo-mechanical analysis of a fractured body with the pressure-dependent thermal contact model; the variations with contact pressure of a the temperature field for the two edges of contact interface, and b the heat flux at the mid-point of the fracture

123

Computational Mechanics

(b) 5

(a)

= 23.86 × 10−4 αa=23.86E-4 = 23.86 × 10−5 αa=23.86E-4 = 23.86 × 10−6 αa=23.86E-4 = 23.86 × 10−7 αa=23.86E-4 = 23.86 × 10−8 αa=23.86E-4 =0 αa=23.86E-4

4

α α α α α α

0

23.86 × 10−4 =a=23.86e-4 =a=23.86e-4 23.86 × 10−5 23.86 × 10−6 =a=23.86e-4 23.86 × 10−7 =a=23.86e-4 23.86 × 10−8 =a=23.86e-4 0 =a=23.86e-4

-0.1

Contact pressure (GPa)

y − displacement (m)

0.1

3

2

2.69 2.690 2.688 2.688 2.686 2.686

1

2.684 2.684 2.682 2.682 2.664 2.664

0 0

1

2

3

4

5

0

1

2

3

2.668 2.688

4

2.672 2.672

5

Load (GPa)

Load (GPa)

Fig. 9 The evolutions of a vertical displacement and b normal contact pressure for various thermal expansion coefficients α at different loading steps

mechanical modeling of a fractured body with frictional contact behavior. This problem was originally proposed by Hirmand et al. [16] to present their augmented Lagrangian contact formulation for non-thermal frictional discontinuity within the X-FEM. The geometry, boundary conditions, and material properties are shown in Fig. 10. The plate is subjected to prescribed displacements in both horizontal and vertical directions as well as a uniform heat flux at the top edge, and it is restrained in two directions and is set to the zero temperature at the bottom edge. The initial temperature value is set to zero for the entire domain. The frictional contact behavior of the fracture surface is studied for two frictional coefficients of μ = 0.1 and 0.4. The thermal contact parameters are chosen as follows; the thermal resistance h c0 = 6.5 × 10−5 N/s◦ C, the Vickers hardness He = 100 N/m2 , and the thermal constant εc = 1. The penalty numbers, which are employed to impose the normal and tangential constraints, are chosen as 1.0×107 MPa. The domain is modeled by a uniform structured FE mesh of 75 × 75 quadrilateral elements. Numerical simulations are performed with different thermal contact conductivity models and the results are compared, including: the pressure-dependent model (PDM) and two simplified conditions, called as the fully conductive model (FCM) and the adiabatic contact model (ACM). In the PDM, the thermal conductivity of the contact surface is related to the contact pressure using a specific constitutive law, as given in relation (15); while in the FCM, the contact surface is modeled with a very high value of thermal conductivity, and in the ACM, it is assumed that there is no heat flux across the contact zone by choosing a very low value of thermal conductivity. It must be emphasized that in the case of an adiabatic contact condition, the thermal boundary condition at the top edge must be changed from the Neu-

123

left Ver

qn

= 0.1 m

right Ver

= 0.01 m

Hor

= 0.05 m

1200 N m2s

fracture interface

hcont E = 1.0 × 104 MPa , ν = 0.3

y

κ = 150 N s o C

x

T =0

Fig. 10 Thermo-mechanical modeling of a fractured body with frictional contact behavior; the geometry, boundary conditions, and material properties

mann type to the Dirichlet type to avoid the singularity in the system of thermal equations; so a prescribed temperature of 10 ◦ C is imposed at the top edge in the adiabatic contact condition. A comparison is finally performed between the results of frictional thermo-mechanical contact analysis using various thermal conductivity models and those obtained from the purely mechanical contact (PMC) analysis with non-thermal condition. In Fig. 11, the evolutions of tangential sliding and normal contact pressure are plotted together with the temperature profile along the fracture surface for two frictional coefficients μ = 0.1 and 0.4. It must be noted that for the frictional coefficient μ = 0.1, the relative tangential displacement reaches the slip limit along the entire contact interface, as

Computational Mechanics

(a)

(b)

0.03

Pressure contact PDM Dependent model ) ( pressure –T dependent Full T contact FCM ( fully conductive model ) No T contact ( adiabatic contact model ) JustACM Mechanical mechanicaletcontact JustPMC Mechanical al. 2015) ( purely(M.Hirmand )

Tangential sliding (m)

Tangential sliding (m)

0.08

0.07

0.06

μ = 0.1 PDM ( pressure – dependent model )

0.05

FCM ( fully conductive model ) Pressure Dependent T contact Full T contact ACM ( adiabatic contact model ) No T contact PMC ( purely mechanical contact ) Just Mechnical PMC et al. 2015 ( Hirmand Just Mechnical (M.Hirmand et)al. 2015)

0.04

0

0.2

0.4

0.6

0.8

0.4

0.6

0.8

1

(d)

PMC ( Hirmand et al. 2015 )

-0.4

0.2

x − coordinate (m)

Pressure Dependent T contact model ) PDM ( pressure – dependent Full T contact FCM ( fully conductive model ) No T contact ACM ( adiabatic contact model ) Just Mechanical PMC contact ( purely mechanical ) Just Mechanical (M.Hirmnd et al. 2015)

-0.2

0.01

0

Temperature ( o C)

Normal contact stress (GPa)

0

μ = 0.4

0 1

x − coordinate (m)

(c)

PMC ( Hirmand et al. 2015 )

0.02

μ = 0.1 & 0.4 -0.6

-0.8

μ = 0.1 & 0.4

Top edge Bottom edge

-1

-1.2 0

0.2

0.4

0.6

0.8

1

x − coordinate (m)

x − coordinate (m)

Fig. 11 The evolutions of tangential sliding and normal contact pressure together with the temperature distribution along the fracture surface for two frictional coefficients μ = 0.1 and 0.4 using different thermal contact conductivity models

can be seen from Fig. 11a; while for the frictional coefficient μ = 0.4, the stick region can be observed between the contacting surfaces where the tangential contact stresses drop below the slip limit, as shown in Fig. 11b. Obviously, the increase of temperature in the domain leads to thermal deformation of the plate in addition to its mechanical deformation that causes additional contact pressure over the contact surface; as a result, the frictional resistance increases and the tangential sliding reduces causing the occurrence of stick zone along the contact surface, as can be seen from Fig. 11a, b. A further comparison between Fig. 11a, b reveals that the contact pressure of the PDM is greater than the FCM model, and accordingly causes the higher value of frictional resistance for the PDM model. This figure clearly illustrates the capability of proposed computational algorithm in capturing the interaction between mechanical and thermal deformations in the frictional thermo-mechanical contact analysis. Moreover, since the normal and tangential contact stresses

are computed in a decoupled manner, the profiles of normal contact pressure obtained from μ = 0.1 are identical to those obtained from μ = 0.4, as shown in Fig. 11c. It can also be seen from Fig. 11d that the results of temperature profile using the pressure-dependent thermal contact conductivity model are identical for both frictional coefficients μ = 0.1 and 0.4 since the thermal conductivity model is dependent on the contact pressure through a non-associated slip rule. Furthermore, it is obvious from Fig. 11d that the thermal contact conductivity decreases along the fracture from the left to right edges of the plate. In Fig. 12, the contours of horizontal displacement are presented on the deformed configurations for frictional coefficients μ = 0.1 and 0.4. A comparison is performed in this figure between different thermal contact conductivity models, i.e. PDM, FCM, and ACM, and the purely mechanical contact (PMC) model with non-thermal condition. In Fig. 13, the contours of temperature distribution are shown

123

Computational Mechanics

Pressure-Dependent Model (PDM)

Fully Conductive Model (FCM)

Adiabatic Contact Model (ACM)

Purely Mechanical Contact (PMC)

(a)

(b) Fig. 12 The contours of horizontal displacement on the deformed configurations; a μ = 0.1, b μ = 0.4; a comparison between different thermal contact conductivity models, i.e. PDM, FCM, and ACM, and the purely mechanical contact (PCM) model with the non-thermal contact condition

Pressure-Dependent Model (PDM)

Fully Conductive Model (FCM)

Adiabatic Contact Model (ACM)

(a)

(b) Fig. 13 The contours of temperature distribution on the deformed configurations; a μ = 0.1, b μ = 0.4 using different thermal contact conductivity models

on the deformed configurations using different thermal contact conductivity models for μ = 0.1 and 0.4. Also presented in Fig. 14 are the contours of heat flux on deformed configurations using the pressure-dependent thermal contact conductivity model for the friction coefficient μ = 0.1, or 0.4. The profile of heat transfer from the top to bottom edges can be clearly observed from the heat flux streamlines, as shown in Fig. 14c. In order to illustrate the performance of the proposed staggered algorithm in comparison with the

123

coupled (monolithic) approach, the evolutions of normalized residual of the displacement and temperature fields are plotted in Fig. 15 at different Newton iterations for friction coefficients μ = 0.1 and 0.4. In this manner, the normalized residual values are defined according to equations (50–b) and (51–b) as ⎛  ⎞   ¯ k+1 −1  ¯ k      dθ ⎜ θ  ⎟  Rtemp = ⎝ ⎠    ω ¯ n   dω  ¯ n+1 

Computational Mechanics

Fig. 14 The contours of heat flux on the deformed configuration using the pressure-dependent thermal contact conductivity model for friction coefficient μ = 0.1, or 0.4; a horizontal direction, b vertical direction, and c the heat flux streamlines

Rdisp

⎛  ⎞   k+1 −1  k      ¯ ¯ d u u ⎜   ⎟  = ⎝  a¯  ⎠  d a¯     n+1 n

(54)

where Rtemp and Rdisp are the normalized residual of the temperature and displacement fields, respectively, and the notation R is the L 2 − norm of vector R. It can be clearly seen that both techniques represent almost similar performance during different Newton iterations. It is noteworthy to emphasize that the proposed staggered approach performs properly the coupling between mechanical and thermal effects.

7.4 Thermo-mechanical modeling of a double-clamped beam with a central crack In this example, a beam with a central crack, which is fixed at both ends, is modeled under thermo-mechanical loading, as shown in Fig. 16. This example was originally proposed by Khoei [43] in purely mechanical condition to determine the active contact zone along the crack interface, where the contact region is a priori unknown. In this study, the doubleclamped beam is simulated to illustrate the performance of the proposed thermo-mechanical computational algorithm in modeling the fracture, which may undergo the opening or closing mode. A uniform structured X-FEM mesh of 3811 bilinear elements is employed, as depicted in Fig. 16. The beam is subjected to the vertical loading of 240 MN at the top edge distributed over a length of 0.4 m at its center. The temperature is set to zero at the left fixed edge and T = 10o C at the right fixed edge. The beam is divided into two parts due to a vertical fracture located at the middle of the beam, and the frictional effects are not considered here since no relative movement occurs along the fracture discontinuity. Hence, the normal contact condition is applied with a penalty parameter of η N = 5×1010 MN/m3 to impose the contact constraint on the crack interface. Since the crack interface can experience

the opening or closing mode, the active contact zone must be determined during the iterative solution along the crack interface. To this end, the entire fracture interface is first discretized by a set of one-dimensional segments; the contact region is then recognized within the Newton–Raphson solution based on a non-positive value of the normal gap distance calculated at the integration points along the crack interface. In Fig. 17, the contours of temperature, horizontal heat flux, horizontal displacement, and stress component σx are presented for different times of loading. Clearly, the results of the final time step of loading shown in Fig. 17e can be compared with those obtained from the steady state analysis, as presented in Fig. 17f. In Fig. 18, the evolutions of the normal contact stress, normal gap displacement, and temperature at two sides of the fracture are plotted along the fracture surface at the final time step of loading. Note that the active contact zone along the fracture surface is obtained within the Newton iterations until the convergence is achieved. It can be seen that the lower part of the fracture surface experiences an opening mode where the normal traction is zero, while the upper part is in contact condition. Obviously, the active contact zone increases in thermo-mechanical contact analysis in comparison to the purely mechanical model, as shown in Fig. 18a, b. Moreover, the contact compression forces are developed in the upper part of the beam, where the closing mode occurs together with the maximum normal contact stresses at the top edge of the beam. The performance of the Heaviside enrichment function within the X-FEM thermo-mechanical contact analysis is clearly obvious from the profiles of temperature at the opening zone along the fracture surface, as shown in Fig. 18c. A comparison is also performed between the results of thermo-mechanical analysis and purely mechanical model in Fig. 19 for the contours of horizontal displacement as well as the stress component σx of the double-clamped beam. Finally, the stability and convergence of the proposed staggered algorithm are presented in Fig. 20 together with the coupled (monolithic) approach. In this figure, the evolutions

123

Computational Mechanics

Fig. 15 The convergence study between the staggered and monolithic approaches in the steady state condition; a, b the normalized residual of the displacement field, c, d the normalized residual of the temperature field

of normalized residual of the displacement and temperature fields are plotted along with the number of iterations at different time steps. The horizontal axis represents the cumulative number of iterations from the beginning to the end of converged solution in the steady-state condition. In this study, the convergence criterion at each time step is set to 1.0E − 4 for the energy norm. Clearly, it can be observed that the residual values of the temperature and displacement fields are less than 1.0E − 8, in which the convergence is obtained with less than three iterations at different time steps of the transient analysis. Obviously, the results of two approaches are comparable during different Newton iterations. It can be seen that the proposed staggered algorithm presents the convergence of the solution with only one iteration at the end of simulation when the steady state condition occurs. It can be highlighted that the proposed staggered approach captures the coupling between mechanical and thermal effects appropriately.

123

7.5 A plate in thermo-mechanical loading condition with multiple cracks The last example is chosen to illustrate the performance of proposed computational algorithm in numerical modeling of a complex fractured geometry for a plate with multiple cracks subjected to thermo-mechanical loading conditions. The plate is modeled employing the X-FEM technique with an unstructured quadrilateral mesh consists of 10,143 elements, as shown in Fig. 21. A linear prescribed displacement with u max = 0.25 m is imposed on the top edge of the plate, and it is fixed in both horizontal and vertical directions over the left and bottom edges. The thermo-mechanical analysis is performed in the steady state condition by applying the fixed temperatures of Th = 100 ◦ C and Tc = 0 at the top and bottom edges of the plate, respectively. The plate contains four inclined cracks with the length of 3m located at different positions and angles, as presented in Table 1. The material properties of the plate are chosen as follows; the elastic mod-

Computational Mechanics Fig. 16 Thermo-mechanical modeling of a double-clamped beam with a central crack; the geometry, boundary conditions, material properties, and X-FEM mesh

0.4 m 240 MN T =0

3m

fracture interface

T = 10o C

c = 100 J kg o C

ρ = 1000 kg m3

κ = 150 N s o C

E = 1× 1010 MPa

β = 23.86 ×106 1 o C

ν = 0.3

0.06 m

10 m

Temperature (°C)

Horizontal flux (N/s.m)

Horizontal displacement (m)

Stress

σ x (Pa)

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 17 The contours of temperature, horizontal heat flux, horizontal displacement, and stress component σx for different times of loading; a t = 1 × 103 s, b t = 3 × 103 s, c t = 5 × 103 s, d t = 1 × 104 s, e t = 2 × 104 s, f the steady state analysis

ulus E = 20 MPa, Poisson ratio v = 0.2, thermal expansion coefficient β = 1.2 × 10−4 1/◦ C, and thermal conductivity κ = 50 N/s ◦ C. The thermal contact properties are assumed as h c0 = 250 N/s ◦ C, He = 45 N/m2 and εc = 1.35. The

tangential sliding of the cracks is set free, and the normal penalty parameter is set to E = 50 GPa. This is an interesting practical example that can be considered as a cross-section of the rock with multiple faults.

123

Computational Mechanics

(c)

y − coordinate (m)

TMC ( thermo − mechanical contact ) PMC ( purely mechanical contact ) PMC ( Khoei 2015 )

TMC ( thermo − mechanical contact ) PMC ( purely mechanical contact ) PMC ( Khoei 2015 )

left side reight side

y − coordinate (m)

(b)

y − coordinate (m)

(a)

Temperature ( o C)

Normal gap (m)

Normal contact stress (MPa)

Fig. 18 The evolutions along the fracture surface of a normal contact stress, b normal gap displacement, and c temperature profile at two sides of fracture at the final time step of loading

Thermo-mechanical contact model

(a)

(b)

-0.07

-8E+08

-6E+08

-0.05

-0.03

-4E+08

-0.01

-2E+08

Purely mechanical contact model

-0.07

0.01

0

2E+08

-8E+08

-6E+08

-0.05

-0.03

-4E+08

-0.01

-2E+08

0.01

0

2E+08

Fig. 19 A comparison between the results of thermo-mechanical analysis and purely mechanical model; a the contours of horizontal displacement, b the contours of stress σx

In order to study the thermal contact behavior in such a complex fractured geometry, it is common to model the thermal contact as an adiabatic condition [17,18,21,56]. In Fig. 22, the contours of horizontal and vertical displacements are shown together with the temperature contours for the thermo-mechanical and adiabatic contact conditions. In this example, the contact pressure is large enough to produce high thermal conductivity at the contact surface, where the heat flux can be transferred perfectly between two crack faces. Moreover, there is no significant effect on the mechanical deformation although different temperature fields can be seen between two thermal contact conditions. Thus, it is essential to model the thermal contact condition appropriately for the temperature-sensitive systems. In Fig. 23, the evolutions of temperature are plotted at the top and bottom faces of crack interface together with the distribution of normal jump along the crack surface for the thermo-mechanical and adiabatic contact models. Obviously, the crack C2 is in complete con-

123

tact condition with zero normal jump of the displacement field (Fig. 23b). In the crack C2 , the continuity of temperature is fully recovered in the thermo-mechanical contact condition, while the temperature of the top and bottom contact faces is different in the adiabatic contact condition. It must be noted that the active contact zones are determined through the Newton–Raphson iterations in both the solid and thermal phases, as different active contact zones can be distinguished for the cracks C1 , C3 and C4 , as shown in Fig. 23. Clearly, it can be observed that the cracks C1 and C4 with larger thermal contact zones have more influence on the displacement field. The importance of thermal contact phenomenon is obvious for this complex fractured geometry. This example clearly illustrates that the proposed computational technique can be used efficiently to model the thermo-mechanical contact problem with less computational efforts for handling of embedded discontinuities.

Computational Mechanics

Fig. 20 The convergence study between the staggered and monolithic approaches for the normalized residual of the displacement and temperature fields at different time steps (the horizontal axis represents the

cumulative number of iterations from the beginning to the end of converged solution in steady-state condition); a t = 102 s, b t = 103 s, c t = 104 s, d t = 105 s

123

Computational Mechanics

max

max

5

Th

4

C3

C4

3

θ4

2 1 0

C2

-1 -2

C1

-3 -4

Tc

-5

-5

-4

-3

-2

-1

0

1

2

3

4

5

(b)

(a)

Fig. 21 A plate with multiple cracks in thermo-mechanical loading conditions; a the geometry, boundary conditions and cracks configuration, b the unstructured quadrilateral mesh consists of 10,143 elements

Table 1 Configuration of the midpoint of crack with the length of 3 m and its direction Crack number 1 2

x-coordinate (m) 1.98 −2

y-coordinate (m)

Angle (degree)

− 2.96

15.0

− 0.96

− 35.0

3

2.51

2.46

55.0

4

− 1.04

2.02

− 75.0

8 Conclusion In the present paper, the extended–FEM technique was employed for modeling the thermo-mechanical contact problem that takes into account the deformable continuum mechanics and the transient heat transfer analysis. The coupled thermo-mechanical analysis was performed by employing the thermal expansion due to thermal effects on the mechanical deformation, and utilizing a nonlinear pressuredependent thermal contact model due to contact forces on the thermal-contact constitute law. To this end, a thermalcontact rule was incorporated within the X-FEM framework as a nonlinear constitutive law that describes the heat flux across the contact zone. The Coulomb frictional law was applied for the mechanical contact problem and a pressuredependent thermal contact model was employed in the weak

123

form of X-FEM method. The Heaviside enrichment function was used for both mechanical and thermal fields to capture the displacement and temperature jumps across the fracture. The governing equations were discretized by the Newmark time splitting method and the final set of nonlinear equations were solved using the Newton–Raphson method in a staggered algorithm. Finally, in order to illustrate the capability and robustness of the proposed computational algorithm several examples were solved numerically, including: the numerical simulation of thermal resistance in a fractured body with constant conductivity at the contact surface, thermo-mechanical modeling of a fractured body with a pressure-dependent thermal contact model, thermomechanical modeling of a fractured body with frictional contact behavior, thermo-mechanical modeling of a doubleclamped beam with a central crack, and modeling of a complex fractured geometry for a plate with multiple cracks. The numerical results were compared with those reported in the literature, or available analytical solutions. In the first example, the performance of thermal contact condition was illustrated in modeling the Kapitza thermal resistance in a fractured body with the constant thermal conductivity at the contact zone. In the second example, the performance of a pressure-dependent thermal contact conductivity model was presented in the thermo-mechanical analysis of a fractured body. It was shown that the assumption of the constant thermal contact conductivity is only acceptable for a spec-

Computational Mechanics

Thermo-mechanical contact Model

Adiabatic Contact Model

(a)

0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2

(b)

(c) Fig. 22 The contours of displacement and temperature distributions for the thermo-mechanical and adiabatic contact models; a horizontal displacement, b vertical displacement, c temperature

ified range of the contact pressure. In the third example, the performance of the thermal-frictional contact behavior was presented in thermo-mechanical modeling of a fractured body using various thermal contact conductivity models, including: the pressure-dependent conductivity model, the fully conductivity model, and the adiabatic contact model. A comparison was performed between the results of various thermal contact conductivity models and those obtained from the purely mechanical contact model with the non-thermal contact condition. The performance of the proposed thermomechanical computational algorithm for a more complex

contact problem was studied through a double-clamped beam in the fourth example to determine the active contact zone along the fracture interface, where the contact region is a priori unknown and may undergo the opening, or closing mode. Finally, the last example was chosen to illustrate the robustness of the proposed computational technique to model the thermo-mechanical contact problem with embedded discontinuities. It was shown that the X-FEM thermo-mechanical contact algorithm proposed here can be used efficiently to model the thermal and mechanical discontinuities along the fracture surface.

123

Computational Mechanics Thermo-mechanical contact Model

Adiabatic Contact Model 0.3

Bottom face Top face Normal jump

40 0.25

0.25

36

0.15 22 0.1

20

Temperature ( C )

0.2 24

Nornal Jump ( cm )

26

Temperature ( C )

0.3

Bottom face Top face Normal jump

0.2

32 28

0.15 24 20

0.1

Normal jump ( cm )

28

44

16 12

16

8

0 1

2

1.5

(a)

2.5

1

1.5

2

2.5

3

3.5

X-coordinate ( m )

X-cordinate ( m )

50

0.3

0.1

42 0

40 38

-0.1

0.2

52

Temperature ( C )

44

Top face Bottom face Normal jump

56

0.2

Normal Jump ( cm )

46

0.3

60

Bottom face Top face Normal jump

48

0.1 48 44

0

40 -0.1 36

36 -0.2

34 32 -2

-1.5

28

-0.3

-1

-3

-2.5

X-coordinate ( m )

-1.5

-1

X-coordinate ( m )

88

0.07

Bottom face Top face Normal jump

0.06

88

0.05

84

0.05

80

0.04

76

0.03

72

0.02

0.01

68

0.01

0

64

0.04 76 0.03 72 0.02

Temperature ( C )

80

0.07

92

Normal jump ( cm )

Bottom face Top face Normal jump

84

-2

0.06

Normal jump

-2.5

-0.2

32

-0.3 -3

(b)

Temperature ( C )

0 0.5

3.5

3

Normal jump ( cm )

0.5

Temperature ( C )

0.05

0.05

18

68

64 1.5

2

2.5

(c)

3

3.5

0 1.5

2

X-coordinate ( m )

1

80

3

3.5

X-coordinate ( m )

Bottom face Top face Normal jump

84

2.5

92

1

Bottom face Top face Normal jump

88

0.8

0.8

0.4

68 64

0.2

80

0.6

76 0.4 72 68

Normal jumo ( cm )

0.6 72

Temperature ( C )

76

Normal jump ( cm )

Temperature ( C )

84

0.2

60 64

56

0

0 60

(d)

-1.4

-1.2

-1

X-coordinate

-0.8

-0.6

-1.4

-1.2

-1

-0.8

-0.6

X-coordinate ( m )

Fig. 23 The evolutions of temperature at the top and bottom faces of crack interface together with the distribution of normal jump along the crack surface for the thermo-mechanical and adiabatic contact models; a crack C1 , b crack C2 , c crack C3 , d crack C4

123

Computational Mechanics

References 1. Khoei AR, Shamloo A, Azami AR (2006) Extended finite element method in plasticity forming of powder compaction with contact friction. Int J Solids Struct 43:5421–5448 2. Khoei AR, Mohammadnejad T (2011) Numerical modeling of multiphase fluid flow in deforming porous media: a comparison between two- and three-phase models for seismic analysis of earth and rockfill dams. Comp Geotech 38:142–166 3. Khoei AR, Hirmand M, Vahab M, Bazargan M (2015) An enriched FEM technique for modeling hydraulically-driven cohesive fracture propagation in impermeable media with frictional natural faults; Numerical and experimental investigations. Int J Numer Methods Eng 104:439–468 4. Khoei AR, Vahab M, Hirmand M (2016) Modeling the interaction between fluid-driven fracture and natural fault using an enrichedFEM technique. Int J Fract 197:1–24 5. Belytschko T, Black T (1999) Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 45:601–620 6. Belytschko T, Moës N, Usui S, Parimi C (2001) Arbitrary discontinuities in finite elements. Int J Numer Methods Eng 50:993–1013 7. Chessa J, Belytschko T (2003) An extended finite element method for two-phase fluids. J Appl Mech 70:10–17 8. Khoei AR, Haghighat E (2011) Extended finite element modeling of deformable porous media with arbitrary interfaces. Appl Math Model 35:5426–5441 9. Dolbow JE, Moës N, Belytschko T (2001) An extended finite element method for modeling crack growth with frictional contact. Comput Methods Appl Mech Eng 190:6825–6846 10. Khoei AR, Nikbakht M (2007) An enriched finite element algorithm for numerical computation of contact friction problems. Int J Mech Sci 49:183–199 11. Liu F, Borja RI (2008) A contact algorithm for frictional crack propagation with the extended finite element method. Int J Numer Meth Eng 76:1489–1512 12. Khoei AR, Vahab M (2014) A numerical contact algorithm in saturated porous media with the extended finite element method. Comput Mech 54:1089–1110 13. Béchet E, Moës N, Wohlmuth B (2009) A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method. Int J Numer Methods Eng 78:931–954 14. Hautefeuille M, Annavarapu C, Dolbow JE (2012) Robust imposition of Dirichlet boundary conditions on embedded surfaces. Int J Numer Methods Eng 90:40–64 15. Annavarapu C, Hautefeuille M, Dolbow JE (2014) A Nitsche stabilized finite element method for frictional sliding on embedded interfaces. Part I: single interface. Comp Methods Appl Mech Eng 268:417–436 16. Hirmand M, Vahab M, Khoei AR (2015) An augmented Lagrangian contact formulation for frictional discontinuities with the extended finite element method. Finite Elem Anal Des 107:28–43 17. Duflot M (2008) The extended finite element method in thermoelastic fracture mechanics. Int J Numer Methods Eng 74:827–847 18. Zamani A, Eslami MR (2010) Implementation of the extended finite element method for dynamic thermoelastic fracture initiation. Int J Solids Struct 47:1392–1404 19. Khoei AR, Moallemi S, Haghighat E (2012) Thermo-hydromechanical modeling of impermeable discontinuity in saturated porous media with X-FEM technique. Eng Fract Mech 96:701– 723 20. Shao Q, Bouhala L, Younes A, Núñez P, Makradi A, Belouettar S (2014) An XFEM model for cracked porous media: effects of fluid flow and heat transfer. Int J Fract 185:155–169

21. Gill P, Davey K (2014) A thermomechanical finite element tool for Leak-before-Break analysis. Int J Numer Methods Eng 98:678–702 22. Yvonnet J, He QC, Toulemonde C (2008) Numerical modelling of the effective conductivities of composites with arbitrarily shaped inclusions and highly conducting interface. Compos Sci Technol 68:2818–2825 23. Yvonnet J, He QC, Zhu QZ, Shao JF (2011) A general and efficient computational procedure for modelling the Kapitza thermal resistance based on XFEM. Comput Mat Sci 50:1220–1224 24. Gu ST, Monteiro E, He QC (2011) Coordinate-free derivation and weak formulation of a general imperfect interface model for thermal conduction in composites. Compos Sci Technol 71:1209–1216 25. Liu JT, Gu ST, Monteiro E, He QC (2014) A versatile interface model for thermal conduction phenomena and its numerical implementation by XFEM. Comput Mech 53:825–843 26. Javili A, McBride A, Steinmann P (2012) Numerical modelling of thermomechanical solids with mechanically energetic (generalised) Kapitza interfaces. Comput Mat Sci 65:542–551 27. Javili A, Kaessmair S, Steinmann P (2014) General imperfect interfaces. Comp Methods Appl Mech Eng 275:76–97 28. Jain A, Tamma KK (2010) Parabolic heat conduction specialized applications involving imperfect contact surfaces: local discontinuous Galerkin finite element method—Part 2. J Therm Stresses 33:344–355 29. Gu ST, Liu JT, He QC (2014) The strong and weak forms of a general imperfect interface model for linear coupled multifield phenomena. Int J Eng Sci 85:31–46 30. Curnier AR, Taylor RL (1982) A thermomechanical formulation and solution of lubricated contacts between deformable solids. J Lubrication Technol 104:109–117 31. Zavarise G, Wriggers P, Stein E, Schrefler BA (1992) Real contact mechanisms and finite element formulation–a coupled thermomechanical approach. Int J Numer Methods Eng 35:767–785 32. Wriggers P, Miehe C (1994) Contact constraints within coupled thermomechanical analysis—a finite element model. Comp Methods Appl Mech Eng 113:301–319 33. Zavarise G, Wriggers P, Schrefler B (1995) On augmented Lagrangian algorithms for thermomechanical contact problems with friction. Int J Numer Methods Eng 38:2929–2949 34. Pantuso D, Bathe KJ, Bouzinov PA (2000) A finite element procedure for the analysis of thermo-mechanical solids in contact. Comput Struct 75:551–573 35. Rieger A, Wriggers P (2004) Adaptive methods for thermomechanical coupled contact problems. Int J Numer Methods Eng 59:871–894 36. Zavarise G, Bacchetto A, Gänser HP (2005) Frictional heating in contact mechanics—a methodology to deal with high temperature gradients. Comput Mech 35:418–429 37. Hansen G (2011) A Jacobian-free Newton Krylov method for mortar-discretized thermo-mechanical contact problems. J Comput Phys 230:6546–6562 38. Hesch C, Franke M, Dittmann M, Temizer I (2016) Hierarchical NURBS and a higher-order phase-field approach to fracture for finite-deformation contact prblems. Comp Methods Appl Mech Eng 301:242–258 39. Belghith S, Mezlini S, Belhadjsalah H, Ligier JL (2013) Thermomechanical modelling of the contact between rough surfaces using homogenisation technique. Mech Res Commun 53:57–62 40. Dittmann M, Franke M, Temizer I, Hesch C (2014) Isogeometric analysis and thermomechanical Mortar contact problems. Comp Methods Appl Mech Eng 274:192–212 41. Murashov MV, Panin SD (2015) Numerical modelling of contact heat transfer problem with work hardened rough surfaces. Int J Heat Mass Trans 90:72–80 42. Grisvard P (2011) Elliptic problems in nonsmooth domains. Vol 69: classics in applied mathematics. SIAM, Philadelphia

123

Computational Mechanics 43. Khoei AR (2015) Extended finite element method: theory and applications. Wiley, New York 44. Le-Quang H, Bonnet G, He QC (2010) Size-dependent Eshelby tensor fields and effective conductivity of composites made of anisotropic phases with highly conducting imperfect interfaces. Phys Rev B 81:064203 45. Benveniste Y, Miloh T (1986) The effective conductivity of composites with imperfect thermal contact at constituent interfaces. Int J Eng Sci 24:1537–1552 46. Lipton R, Vernescu B (1996) Composites with imperfect interface. Proc R Soc Lond A Math Phys Eng Sci. https://doi.org/10.1098/ rspa.1996.0018 47. Hashin Z (2001) Thin interphase/imperfect interface in conduction. J Appl Phys 89:2261–2267 48. Benveniste Y, Miloh T (1999) Neutral inhomogeneities in conduction phenomena. J Mech Phys Solids 47:1873–1892 49. Madhusudana CV, Fletcher LS (1986) Contact heat transfer—the last decade. AIAA J 24:510–523 50. Grosch KA (1963) The relation between the friction and viscoelastic properties of rubber. Proc R Soc Lond A 274:21–39 51. Mase CW, Smith L (1987) Effect of frictional heating on the thermal, hydrologic, and mechanical response of a fault. J Geophys Res 92:6249–6272

123

52. Braun OM, Steenwyk B, Warhadpande A, Persson BNJ (2016) On the dependency of friction on load: theory and experiment. Europhys Lett 113:56002 53. Farhat C, Park KC, Yves DP (1991) An unconditionally stable staggered algorithm for transient finite element analysis of coupled thermoelastic problems. Comput Methods Appl Mech Eng 85:349– 365 54. Turska E, Schrefler BA (1993) On convergence conditions of partitioned solution procedures for consolidation problems. Comput Methods Appl Mech Eng 106:51–63 55. Danowski C, Gravemeier V, Yoshihara L, Wall WA (2013) A monolithic computational approach to thermo-structure interaction. Int J Numer Methods Eng 95:1053–1078 56. Nguyen MN, Bui TQ, Nguyen NT, Truong TT, Lich LV (2017) Simulation of dynamic and static thermoelastic fracture problems by extended nodal gradient finite elements. Int J Mech Sci 134:370– 386 57. Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T (1999) Computational geomechanics with special reference to earthquake engineering. Wiley, New York