- Email: [email protected]
Symmetrical node-to-node formulation for thermal contact problems between non-conforming or non-matching meshes
Finite Elements in Analysis and Design 162 (2019) 13–18
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Symmetrical node-to-node formulation for thermal contact problems between non-conforming or non-matching meshes Eric Feulvarch a, ∗ , Jean-Christophe Roux a , Marion Geuffrard b a b
Univ. Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 Rue Jean Parot, 42023 Saint-Etienne Cedex 2, France ESI Group, Le Recamier, 70 Rue Robert, 69458 Lyon Cedex 06, France
A R T I C L E
I N F O
Keywords: Thermal contact Finite element Non-conforming meshes Non-matching meshes
A B S T R A C T
The aim of this paper is to propose a symmetrical finite element formulation for the numerical treatment of thermal contact problems without any constraint on the shape, size and order of approximation of the meshes in contact. The new formulation implicitly ensures the heat balance through the contact and it is well adapted for the use of non-linear symmetric solvers minimizing the matrix storage. Examples are presented to show the relevance of the formulation developed.
1. Introduction The numerical simulation of heat transfer often leads to simulating contact problems by using a thermal contact resistance or a contact heat transfer coefficient. This is particularly true in the field of manufacturing processes. One can cite the electrode/sheet contact in spot welding, the sheet/backing plate contact in friction stir welding or the insert/tool holder contact in machining. A simple FEM approach consists in meshing the contact surfaces in the same way as proposed by Wriggers and Zavaris [1] and Johansson and Klarbring [2] who consider that the heat is transferred between nodes that are face to face. In the same way, Robin et al. [3] proposed an element-to-element formulation similar to the one developed by Pantuso et al. [4] where the heat is exchanged between elements that are face to face. It has been successfully used by Raoelison et al. [5] for simulating the resistance spot welding process. Unfortunately, the meshes of components are often build separately and thus, are not identical. Moreover, the possible mechanical movements can change the geometrical configuration (ex: hot stamping [6] or spark plasma sintering [7]) and the meshes that were initially identical can become very different in terms of element size and integration point connectivity. In this case, Agelet de Saracibar [8], Laursen [9] and Oancea et al. [10] have proposed different coupled thermo-mechanical formulations. As detailed by Dittmann et al. [11], Hueber et al. [12], Molinari et al. [13], Seitz et al. [14] and Temizer [15], the numerical formulations developed are
not natively symmetric for the thermal aspects in the case of nonmatching or non-conforming meshes. The formulations proposed by Stromber [16] and Wriggers and Miehe [17] are only available in 2D. A FEM approach has been proposed by Feulvarch et al. [18] in order to avoid the mesh constraints. This is an integration point-toelement formulation similar to the one developed by Xing and Makinouchi [19] that consists to estimate the heat flux at an integration point by using the temperature of its projection on the opposite element. In the same way, it is possible to build a node-to-element approach as it is done for contact mechanics in some commercial finite element codes [20,21]. Unfortunately, such approaches leads again to a nonsymmetrical finite element formulation that can not be used with symmetric solvers minimizing the matrix storage. Another drawback of the integration point-to-element formulation is the unrealistic results that can be obtained for high values of the contact heat transfer coefficient. This phenomenon is due to the multiple connectivity of nodes across the contact surface induced by the heat exchanged calculation at the integration points. Trujillo and Pappoff [22] have proposed another general FEM formulation that is symmetric but this requires special attention for the numerical integration. The authors proposed examples with one surface having always a finer mesh than the opposite surface. The purpose of this paper is to propose a new node-to-node formulation for the numerical treatment of thermal contact problems without any constraint on:
∗ Corresponding author. E-mail address: [email protected] (E. Feulvarch). https://doi.org/10.1016/j.finel.2019.05.003 Received 8 November 2018; Received in revised form 5 April 2019; Accepted 8 May 2019 Available online 20 May 2019 0168-874X/© 2019 Elsevier B.V. All rights reserved.
E. Feulvarch et al.
• • • •
Finite Elements in Analysis and Design 162 (2019) 13–18
the choice of the solver that can be symmetric, the value of the contact heat transfer coefficient, the shapes and the sizes of the meshes that are face to face, the orders of approximation for each mesh in contact.
Contrary to the work of Jin et al. [23] dedicated to contact mechanics, there is no adding of node in the symmetrical formulation that is described in the first part of the paper. It does not need the control of the spatial discretization as proposed by Rieger and Wriggers [24]. It considers that the heat is exchanged between two nodes that are on both sides of the interface. The balance of heat exchanged is ensured by considering that heat is passing through a unique intermediate mean surface. This interface defines the size of the integration domain for all nodes belonging to both contact surfaces. Examples are presented in the last part to show the relevance of the formulation developed.
Fig. 1. First step for building the nodal connectivities across the contact surface: the nodes of the surface A are connected with the nearest nodes on the surface B.
2. Numerical formulation 2.1. The thermal contact problem Considering the Fourier conduction law, the heat equation writes: {
𝜌 C 𝜃̇ 𝜆 grad 𝜃 · n
in Ω = div (𝜆 grad 𝜃) = qext + hext (𝜃 ext − 𝜃) on 𝜕Ω
(1)
where 𝜌 is the material volumetric mass; C is the specific heat; 𝜃 is the temperature; 𝜆 is the thermal conductivity; n is the unit outward normal vector to the boundary ∂Ω of the domain Ω; qext is a prescribed input heat flux; hext is a heat transfer coefficient for modeling the heat exchanged between the component studied and the external media supposed to be at the temperature 𝜃 ext . As is well-known, the general boundary conditions (1)2 encompass both the cases of a prescribed flux (for hext = 0) and a prescribed value of 𝜃 (for hext → + ∞). The weak formulation of the heat equation is classically obtained by multiplying Eq. (1)1 by weighing function 𝜃 ∗ and integrating over the domain Ω [25]. Integrating by parts and accounting for the boundary condition on the boundary ∂Ω, one thus obtains the following weak formulation: ∫Ω
grad𝜃 ∗ 𝜆 grad𝜃 dv +
=
(
∫𝜕Ω
∫Ω
Fig. 2. Second step for building the nodal connectivities across the contact surface: local association of ‘free’ nodes on the surface B with their nearest nodes on the surface A . B
for the first surface A and the second surface B with nA and nB nodes respectively. The connectivity algorithm is divided into 2 steps. For each node A B a ∈ , we first look for the nearest node b ∈ and we connect them by creating an edge (a, b) and thus a bipartite graph1 between A B nodes of and . One can note that, at this stage, all the nodes B in belonging to the surface B are not necessarily connected to the nodes of the surface A (see Fig. 1). To achieve the complete connectivity across the contact surface for all nodes belonging to the contact, a second loop is carried out for the B ‘free’ nodes in on the surface B , applying the following rule 1. Rule 1. (see Fig. 2). For each ‘free’ node b ∈ NB that has not been connected to one node A A in in the first loop, we look for its nearest node a ∈ . If a is already B connected (after the first loop) to another node c ∈ , and if c is conA nected to at least one more node d in different from node a, then edge (a, c) is removed (i.e. a and c are disconnected). As a consequence of this rule, there exists no path of length greater than two in the bipartite graph. This limits the connections between the far nodes of each surfaces. The objective of this approach is to build only ‘local’ connectivities as shown in Fig. 3 for minimizing the numerical stiffness. In the case of a perfect contact condition, the nodes a and b in Fig. 2 will have the same temperature that could be different than the one of the nodes c and d. Another interesting and useful consequence of this approach is that it generates nC connected components (i.e. nC subsets) Ck (1 ⩽ k ⩽ nC ) of ‘locally’ connected nodes (see Fig. 3). This allows us to define an indicator function relative to each connected component k. This function indicates that two nodes are connected inside a same component:
𝜃 ∗ 𝜌 C 𝜃̇ dv )
𝜃 ∗ qext + hext (𝜃 ext − 𝜃) ds
(2)
In the case of a contact between two component A and B, the righthand side of Eq. (2) leads to the following expressions of the boundary conditions for finite element nodes a and b belonging to the contact surfaces A and B respectively: (
∫𝜕ΩA ⧵ A
)
Ψa qext + hext (𝜃 ext − 𝜃) ds −
Ψa hct (𝜃 A − 𝜃 B ) ds
(3)
Ψb hct (𝜃 A − 𝜃 B ) ds
(4)
∫ A
and (
∫𝜕ΩB ⧵ B
)
Ψb qext + hext (𝜃 ext − 𝜃) ds +
∫ B
where Ψa and Ψb denote respectively the shape functions associated to the a and b; hct is the contact heat transfer coefficient. Thus, the aim of this work is to develop a new numerical method for computing the integrals ∫ A Ψa h (𝜃 B − 𝜃 A ) ds and ∫ B Ψb h (𝜃 A − 𝜃 B ) ds by ensuring the heat balance at the interface. 2.2. Nodal connectivities across the contact surface
1
A bipartite graph (or biograph) is a graph whose vertices can be divided into A B two disjointed and independent sets and such that each edge connects A B a vertex in to one in . Equivalently, a bipartite graph is a graph that does not contain any odd-length cycle.
At first, the nodes of each contact surface are selected by eliminating the nodes that have no orthogonal projection on the opposite surface A and the remaining nodes form two numbered sets of nodes and 14
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 3. Connected components of the final bipartite graph of connectivities between nodes across the contact surface. A
Fig. 4. Sketch of the thermal problem studied.
S∕SA is a global corrective term that ensures the correct computation of the global size of the heat exchanged surface S:
B
Definition 1. For any two nodes (a, b) ∈ × and any connected component Ck , A
a∈
A
(5)
3. Applications 3.1. Perfect contact condition
∑ A
∫ A
B
∫ B
Ψa ds
(6) Tests are performed on the problem defined in Fig. 4. The value of the thermal conductivity is taken equal to 40W∕(m.◦ C). The heat transfer coefficient is 106 W∕(m2 .◦ C) on the left side to prescribe a temperature 𝜃 equal to the external temperature 𝜃 ext = 20 ◦ C. The heat transfer coefficient equals 200W∕(m2 .◦ C) on the right side and the external temperature is of 500 ◦ C. The contact surfaces A and B are chosen to be not perpendicular to the heat transfer direction to highlight the efficiency of the formulation proposed. The simulations have been performed in a static way by using the computer code Sysweld® [21]. At first, we study a perfect thermal contact problem by considering a very high contact heat transfer coefficient hct of 109 W∕(m2 .◦ C) in Eq. (8). This situation corresponds to the most difficult one to compute because it corresponds to the gluing of non-conforming meshes. As a reference, the exact solution is plotted on a global mesh made of P1 elements on the whole geometry in Fig. 5. By considering the contact between two meshes made of P1 elements which are not face to face, a first test is performed using the
and ∑
Ψb ds
(7)
All integrals are computed by means of the classical Gaussian quadrature. 2.4. Nodal calculation of heat exchanged The new formulation consists of evaluating the heat exchanged A B between two connected nodes (a, b) ∈ × (corresponding to A B A ∫ A Ψa hct (𝜃 − 𝜃 ) ds in Eq. (3) and ∫ B Ψb hct (𝜃 − 𝜃 B ) ds in Eq. (4)) by a unique quantity Q(a,b) defined as follows: Q(a,b) = hct (𝜃a − 𝜃b ) S(a,b)
(8)
where 𝜃 a and 𝜃 b denote the temperatures of the two nodes a and b belonging respectively to the surfaces A and B ; hct is the heat exchange coefficient; S(a,b) is the nodal integration surface defined by: ⎛ ⎞ S ⎜ 1ICk (a, b) ∫ B Ψb ds ⎟ S(a,b) = A ⎜ ∑ Ψ ds (9) S ⎜ 1ICk (a, c) ∫ B Ψc ds ⎟⎟ ∫ A a ⎝ c ∈ B ⎠ ∑ The term 1ICk (a, b) ∫ B Ψb ds ∕ c∈ B 1ICk (a, c) ∫ B Ψc ds is a partition coefficient different than 1 when the node a is connected to more than one node b. It is such that: A
(11)
can also be estimated by exchanging the two surfaces A and B but its value will be modified. As for the connectivity procedure (see Remark 1), the finite element formulation proposed is natively not commutative but it always leads to a symmetric equations system.
with
b ∈
a∈
In Eq. (9), the surface A is taken as the integration support. S(a,b)
1 (S + S B ) 2 A
SB =
B
S Ψ ds = S S A ∫ A a
received by B . Therefore, the finite element formulation proposed is symmetric.
The size of the surface of heat exchanged is evaluated by considering the mean values of the parts of each surface A and B belonging to the contact. The mean surface is computed as follows:
a∈
b ∈
∑
Remark 2. The balance of heat exchanged is implicitly ensured because Q(a,b) estimates simultaneously the heat lost locally by A and the one
2.3. Surface of heat exchanged
SA =
A
S(a,b) =
In order to reduce the computing time, S(a,b) can be evaluated only one time once the mesh is build.
Remark 1. The connectivity procedure proposed is not necessary commutative because of the rule 1. Indeed, the subsets Ck are not necessary the same if the surfaces SA and SB are exchanged.
S=
∑
∑
B
1ICk : × → {0, 1} { A B 1 if {a, b} ∈ Ck and (a, b) ∈ × , 1ICk (a, b) ≔ 0 otherwise.
∀a ∈ ,
∑ b ∈
∑ B
1ICk (a, b) ∫ B Ψb ds
c ∈
B
1ICk (a, c) ∫ B Ψc ds
=1
(10)
Fig. 5. 2D reference temperature distribution (◦ C). 15
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 8. 3D temperature distribution (◦ C) between 2 meshes of 6 mm thickness made of P1 (left) and P2 (right) finite elements.
Fig. 6. 2D temperature distribution (◦ C) with P1 elements.
new formulation developed in this work and the technique proposed by Feulvarch et al. [18]. The results are satisfactory for the new technique whereas the temperature distribution is unrealistic for the second method (see Fig. 6). The contact seems to be too stiff because the temperature is constant on all the contact surfaces. Different 2D tests have been carried out on Fig. 7 by using higher order P2 and Q2 elements for the mesh of the part located on the right. One can note that the results are very satisfactory despite the high difference in the mesh sizes. Fig. 8 shows the 3D temperature distribution in the case a contact between a mesh made of P1 elements and a second one made of larger P2 elements. This calculation shows once again the relevance of the formulation proposed in this paper for 3D applications.
Fig. 9. Sketch of the thermal problem studied.
3.2. Convergence test This application is similar to the one of the previous application with one added hole as shown in Fig. 9 for having a non-linear spatial evolution of the temperature. We propose to test the convergence of the numerical solution 𝜃 h obtained with P1 finite elements to a reference solution 𝜃 plotted in Fig. 10. The field 𝜃 was built by means of a computation on a 2D mesh made of finite elements of type P1 ◦ (h = 2.10−2 mm) without considering contact conditions. In Fig. 11, the order of convergence 𝛼 is identified by means of the slope of the linear regression on the log-log plot of ‖𝜃 − 𝜃 h ‖L2 (Ω) (≈C hα ). Two first tests have been performed by considering a homogeneous element size on the whole mesh. For the first one denoted T1 , 1 and 2 correspond respectively to A and B whereas for the second test T2 , A and B are exchanged. Then, two other tests have been carried out by considering a fix element size on the left part of the mesh equal to the one used to build the reference solution. The element size h only evolves on the right part of the mesh. As above, 1
Fig. 7. 2D temperature distribution (◦ C) for different types of meshes.
Fig. 10. 2D reference temperature distribution (◦ C). 16
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 11. Convergence plot of the numerical solution 𝜃 h .
and 2 correspond respectively to A and B for the test T3 whereas they are exchanged for the test T4 . One can note that a first order accuracy is obtained for the tests T1 and T2 corresponding to a homogeneous evolution of the element size h. Moreover, the permutation of the contact surfaces has no significant effect on the quality of the solution and the convergence rate. For the tests T3 and T4 , it is difficult to analyse the value of the convergence
Fig. 13. 3D temperature distribution (o C) between 2 meshes of 6 mm thickness made of P1 (left) and P2 (right) finite elements.
rate because there is not a homogeneous evolution of the element size. The interest of these tests is to show that the permutation of the contact surfaces does not affect the quality and the convergence of the solution in the case of non-matching meshes. 3.3. Contact with thermal contact resistance In this section, we consider a contact heat transfer coefficient hct of 103 W∕(m2 .◦ C) in Eq. (8). Then, the temperature field is not continuous from one side to the other side of the contact. The 2D temperature distribution is plotted in Fig. 12 for the contact between a mesh made of P1 elements on the left and a mesh made of P2 or Q2 elements on the right. The corresponding distribution in 3D is shown in Fig. 13. Once again, the results are very satisfactory compared to the reference solution in Fig. 12. 4. Conclusion The aim of this work was to propose a new formulation for the numerical treatment of thermal contact problems without any constraint on the discretization of the contact surfaces and on the order of the meshes in contact. The results are very satisfactory for all the examples and they show the relevance of the formulation developed for thermal contact problems. References [1] P. Wriggers, G. Zavaris, Thermomechanical contact - a rigorous but simple numerical approach, Comput. Struct. 46 (1993) 47–53. [2] L. Johansson, A. Klarbring, Thermoelastic frictional contact problems: modelling, finite element approximation and numerical realization, Comput. Methods Appl. Mech. Eng. 105 (1993) 181–210. [3] V. Robin, A. Sanchez, T. Dupuy, J. Soigneux, J.M. Bergheau, Numerical simulation of spot welding with special attention to contact conditions, Graz-Seggau, Austria, October 2001, in: H. Cerjak (Ed.), Mathematical Modelling of Weld Phenomena 6, The Institute of Materials, London, 2002. [4] D. Pantuso, K.-J. Bathe, P.A. Bouzinov, A finite element procedure for the analysis of thermo-mechanical solids in contact, Comput. Struct. 75 (6) (2000) 551–573. [5] R.N. Raoelison, A. Fuentes, C. Pouvreau, Ph Rogeon, P. Carre, F. Dechalotte, Modeling and numerical simulation of the resistance spot welding of zinc coated steel sheets using rounded tip electrode: analysis of required conditions, Appl. Math. Model. 38 (2014) 2505–2521. [6] Mohsen Loh-Mousavi, Mehdi AhmadiRad, Tomoyoshi Maeno, Denis J. Politis, LiLiang Wang, Coupled thermal-electrical finite element analysis of electrical resistance heating in hot stamping of ultra-high strength steel tubes, Procedia Manuf. 15 (2018) 1047–1054.
Fig. 12. 2D temperature distribution (◦ C) for different types of meshes. 17
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18 [16] N. Stromber, Finite element treatment of two-dimensional thermoelastic wear problems, Comput. Methods Appl. Mech. Eng. 177 (3) (1999) 441–455. [17] P. Wriggers, C. Miehe, Contact constraints within coupled thermomechanical analysis - a finite element model, Comput. Methods Appl. Mech. Eng. 113 (1994) 301–319. [18] E. Feulvarch, V. Robin, J.M. Bergheau, Resistance spot welding simulation: a general finite element formulation of electrothermal contact conditions, J. Mater. Process. Technol. 153154 (2004) 436–441. [19] H. Xing, A. Makinouchi, Three dimensional finite element modeling of thermomechanical frictional contact between finite deformation bodies using R-minimum strategy, Comput. Methods Appl. Mech. Eng. 191 (3738) (2002) 4193–4214. [20] ANSYS, Users Manual, 2018. [21] ESI Group, Users Manual, 2018. [22] David M. Trujillo, Chris G. Pappoff, A general thermal contact resistance nite element, Finite Elem. Anal. Des. 38 (2002) 263–276. [23] Seungmin Jin, Dongwoo Sohn, Im Seyoung, Node-to-node scheme for three-dimensional contact mechanics using polyhedral type variable-node elements, Comput. Methods Appl. Mech. Eng. 304 (2016) 217–242. [24] A. Rieger, P. Wriggers, Adaptive methods for thermomechanical coupled contact problems, Int. J. Numer. Methods Eng. 59 (2004) 871–894. [25] E. Feulvarch, J.M. Bergheau, J.B. Leblond, An implicit finite element algorithm for the simulation of diffusion with phase changes in solids, Int. J. Numer. Methods Eng. 78 (2009) 1492–1512.
[7] C. Maniere, A. Pavia, L. Durand, G. Chevallier, K. Afanga, C. Estournes, Finite-element modeling of the electro-thermal contacts in the spark plasma sintering process, J. Eur. Ceram. Soc. 36 (2016) 741–748. [8] C. Agelet de Saracibar, Numerical analysis of coupled thermomechanical frictional contact problems. computational model and applications, Arch. Comput. Methods Eng. 5 (3) (1998) 243–301. [9] T.A. Laursen, On the development of thermodynamically consistent algorithms for thermomechanical frictional contact, Comput. Methods Appl. Mech. Eng. 177 (3) (1999) 273–287. [10] V.G. Oancea, T.A. Laursen, A finite element formulation of thermomechanical rate-dependent frictional sliding, Int. J. Numer. Methods Eng. 40 (23) (1997) 4275–4311. [11] M. Dittmann, M. Franke, I. Temizer, C. Hesch, Isogeometric analysis and thermomechanical mortar contact problems, Comput. Methods Appl. Mech. Eng. 274 (2014) 192–212. [12] S. Hueber, B.I. Wohlmuth, Thermo-mechanical contact problems on non-matching meshes, Comput. Methods Appl. Mech. Eng. 198 (15) (2009) 1338–1350. [13] J. Molinari, M. Ortiz, R. Radovitzky, E. Repetto, Finite-element modeling of dry sliding wear in metals, Eng. Comput. 18 (2001) 592–610. [14] A. Seitz, W.A. Wall, A. Popp, A computational approach for thermo-elasto-plastic frictional contact based on a monolithic formulation using non-smooth nonlinear complementarity functions, Adv. Model. Simul. Eng. Sci 5 (1) (2018) 5. [15] I. Temizer, Sliding friction across the scales: thermomechanical interactions and dissipation partitioning, J. Mech. Phys. Solids 89 (2016) 126–148.
18
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Symmetrical node-to-node formulation for thermal contact problems between non-conforming or non-matching meshes Eric Feulvarch a, ∗ , Jean-Christophe Roux a , Marion Geuffrard b a b
Univ. Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 Rue Jean Parot, 42023 Saint-Etienne Cedex 2, France ESI Group, Le Recamier, 70 Rue Robert, 69458 Lyon Cedex 06, France
A R T I C L E
I N F O
Keywords: Thermal contact Finite element Non-conforming meshes Non-matching meshes
A B S T R A C T
The aim of this paper is to propose a symmetrical finite element formulation for the numerical treatment of thermal contact problems without any constraint on the shape, size and order of approximation of the meshes in contact. The new formulation implicitly ensures the heat balance through the contact and it is well adapted for the use of non-linear symmetric solvers minimizing the matrix storage. Examples are presented to show the relevance of the formulation developed.
1. Introduction The numerical simulation of heat transfer often leads to simulating contact problems by using a thermal contact resistance or a contact heat transfer coefficient. This is particularly true in the field of manufacturing processes. One can cite the electrode/sheet contact in spot welding, the sheet/backing plate contact in friction stir welding or the insert/tool holder contact in machining. A simple FEM approach consists in meshing the contact surfaces in the same way as proposed by Wriggers and Zavaris [1] and Johansson and Klarbring [2] who consider that the heat is transferred between nodes that are face to face. In the same way, Robin et al. [3] proposed an element-to-element formulation similar to the one developed by Pantuso et al. [4] where the heat is exchanged between elements that are face to face. It has been successfully used by Raoelison et al. [5] for simulating the resistance spot welding process. Unfortunately, the meshes of components are often build separately and thus, are not identical. Moreover, the possible mechanical movements can change the geometrical configuration (ex: hot stamping [6] or spark plasma sintering [7]) and the meshes that were initially identical can become very different in terms of element size and integration point connectivity. In this case, Agelet de Saracibar [8], Laursen [9] and Oancea et al. [10] have proposed different coupled thermo-mechanical formulations. As detailed by Dittmann et al. [11], Hueber et al. [12], Molinari et al. [13], Seitz et al. [14] and Temizer [15], the numerical formulations developed are
not natively symmetric for the thermal aspects in the case of nonmatching or non-conforming meshes. The formulations proposed by Stromber [16] and Wriggers and Miehe [17] are only available in 2D. A FEM approach has been proposed by Feulvarch et al. [18] in order to avoid the mesh constraints. This is an integration point-toelement formulation similar to the one developed by Xing and Makinouchi [19] that consists to estimate the heat flux at an integration point by using the temperature of its projection on the opposite element. In the same way, it is possible to build a node-to-element approach as it is done for contact mechanics in some commercial finite element codes [20,21]. Unfortunately, such approaches leads again to a nonsymmetrical finite element formulation that can not be used with symmetric solvers minimizing the matrix storage. Another drawback of the integration point-to-element formulation is the unrealistic results that can be obtained for high values of the contact heat transfer coefficient. This phenomenon is due to the multiple connectivity of nodes across the contact surface induced by the heat exchanged calculation at the integration points. Trujillo and Pappoff [22] have proposed another general FEM formulation that is symmetric but this requires special attention for the numerical integration. The authors proposed examples with one surface having always a finer mesh than the opposite surface. The purpose of this paper is to propose a new node-to-node formulation for the numerical treatment of thermal contact problems without any constraint on:
∗ Corresponding author. E-mail address: [email protected] (E. Feulvarch). https://doi.org/10.1016/j.finel.2019.05.003 Received 8 November 2018; Received in revised form 5 April 2019; Accepted 8 May 2019 Available online 20 May 2019 0168-874X/© 2019 Elsevier B.V. All rights reserved.
E. Feulvarch et al.
• • • •
Finite Elements in Analysis and Design 162 (2019) 13–18
the choice of the solver that can be symmetric, the value of the contact heat transfer coefficient, the shapes and the sizes of the meshes that are face to face, the orders of approximation for each mesh in contact.
Contrary to the work of Jin et al. [23] dedicated to contact mechanics, there is no adding of node in the symmetrical formulation that is described in the first part of the paper. It does not need the control of the spatial discretization as proposed by Rieger and Wriggers [24]. It considers that the heat is exchanged between two nodes that are on both sides of the interface. The balance of heat exchanged is ensured by considering that heat is passing through a unique intermediate mean surface. This interface defines the size of the integration domain for all nodes belonging to both contact surfaces. Examples are presented in the last part to show the relevance of the formulation developed.
Fig. 1. First step for building the nodal connectivities across the contact surface: the nodes of the surface A are connected with the nearest nodes on the surface B.
2. Numerical formulation 2.1. The thermal contact problem Considering the Fourier conduction law, the heat equation writes: {
𝜌 C 𝜃̇ 𝜆 grad 𝜃 · n
in Ω = div (𝜆 grad 𝜃) = qext + hext (𝜃 ext − 𝜃) on 𝜕Ω
(1)
where 𝜌 is the material volumetric mass; C is the specific heat; 𝜃 is the temperature; 𝜆 is the thermal conductivity; n is the unit outward normal vector to the boundary ∂Ω of the domain Ω; qext is a prescribed input heat flux; hext is a heat transfer coefficient for modeling the heat exchanged between the component studied and the external media supposed to be at the temperature 𝜃 ext . As is well-known, the general boundary conditions (1)2 encompass both the cases of a prescribed flux (for hext = 0) and a prescribed value of 𝜃 (for hext → + ∞). The weak formulation of the heat equation is classically obtained by multiplying Eq. (1)1 by weighing function 𝜃 ∗ and integrating over the domain Ω [25]. Integrating by parts and accounting for the boundary condition on the boundary ∂Ω, one thus obtains the following weak formulation: ∫Ω
grad𝜃 ∗ 𝜆 grad𝜃 dv +
=
(
∫𝜕Ω
∫Ω
Fig. 2. Second step for building the nodal connectivities across the contact surface: local association of ‘free’ nodes on the surface B with their nearest nodes on the surface A . B
for the first surface A and the second surface B with nA and nB nodes respectively. The connectivity algorithm is divided into 2 steps. For each node A B a ∈ , we first look for the nearest node b ∈ and we connect them by creating an edge (a, b) and thus a bipartite graph1 between A B nodes of and . One can note that, at this stage, all the nodes B in belonging to the surface B are not necessarily connected to the nodes of the surface A (see Fig. 1). To achieve the complete connectivity across the contact surface for all nodes belonging to the contact, a second loop is carried out for the B ‘free’ nodes in on the surface B , applying the following rule 1. Rule 1. (see Fig. 2). For each ‘free’ node b ∈ NB that has not been connected to one node A A in in the first loop, we look for its nearest node a ∈ . If a is already B connected (after the first loop) to another node c ∈ , and if c is conA nected to at least one more node d in different from node a, then edge (a, c) is removed (i.e. a and c are disconnected). As a consequence of this rule, there exists no path of length greater than two in the bipartite graph. This limits the connections between the far nodes of each surfaces. The objective of this approach is to build only ‘local’ connectivities as shown in Fig. 3 for minimizing the numerical stiffness. In the case of a perfect contact condition, the nodes a and b in Fig. 2 will have the same temperature that could be different than the one of the nodes c and d. Another interesting and useful consequence of this approach is that it generates nC connected components (i.e. nC subsets) Ck (1 ⩽ k ⩽ nC ) of ‘locally’ connected nodes (see Fig. 3). This allows us to define an indicator function relative to each connected component k. This function indicates that two nodes are connected inside a same component:
𝜃 ∗ 𝜌 C 𝜃̇ dv )
𝜃 ∗ qext + hext (𝜃 ext − 𝜃) ds
(2)
In the case of a contact between two component A and B, the righthand side of Eq. (2) leads to the following expressions of the boundary conditions for finite element nodes a and b belonging to the contact surfaces A and B respectively: (
∫𝜕ΩA ⧵ A
)
Ψa qext + hext (𝜃 ext − 𝜃) ds −
Ψa hct (𝜃 A − 𝜃 B ) ds
(3)
Ψb hct (𝜃 A − 𝜃 B ) ds
(4)
∫ A
and (
∫𝜕ΩB ⧵ B
)
Ψb qext + hext (𝜃 ext − 𝜃) ds +
∫ B
where Ψa and Ψb denote respectively the shape functions associated to the a and b; hct is the contact heat transfer coefficient. Thus, the aim of this work is to develop a new numerical method for computing the integrals ∫ A Ψa h (𝜃 B − 𝜃 A ) ds and ∫ B Ψb h (𝜃 A − 𝜃 B ) ds by ensuring the heat balance at the interface. 2.2. Nodal connectivities across the contact surface
1
A bipartite graph (or biograph) is a graph whose vertices can be divided into A B two disjointed and independent sets and such that each edge connects A B a vertex in to one in . Equivalently, a bipartite graph is a graph that does not contain any odd-length cycle.
At first, the nodes of each contact surface are selected by eliminating the nodes that have no orthogonal projection on the opposite surface A and the remaining nodes form two numbered sets of nodes and 14
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 3. Connected components of the final bipartite graph of connectivities between nodes across the contact surface. A
Fig. 4. Sketch of the thermal problem studied.
S∕SA is a global corrective term that ensures the correct computation of the global size of the heat exchanged surface S:
B
Definition 1. For any two nodes (a, b) ∈ × and any connected component Ck , A
a∈
A
(5)
3. Applications 3.1. Perfect contact condition
∑ A
∫ A
B
∫ B
Ψa ds
(6) Tests are performed on the problem defined in Fig. 4. The value of the thermal conductivity is taken equal to 40W∕(m.◦ C). The heat transfer coefficient is 106 W∕(m2 .◦ C) on the left side to prescribe a temperature 𝜃 equal to the external temperature 𝜃 ext = 20 ◦ C. The heat transfer coefficient equals 200W∕(m2 .◦ C) on the right side and the external temperature is of 500 ◦ C. The contact surfaces A and B are chosen to be not perpendicular to the heat transfer direction to highlight the efficiency of the formulation proposed. The simulations have been performed in a static way by using the computer code Sysweld® [21]. At first, we study a perfect thermal contact problem by considering a very high contact heat transfer coefficient hct of 109 W∕(m2 .◦ C) in Eq. (8). This situation corresponds to the most difficult one to compute because it corresponds to the gluing of non-conforming meshes. As a reference, the exact solution is plotted on a global mesh made of P1 elements on the whole geometry in Fig. 5. By considering the contact between two meshes made of P1 elements which are not face to face, a first test is performed using the
and ∑
Ψb ds
(7)
All integrals are computed by means of the classical Gaussian quadrature. 2.4. Nodal calculation of heat exchanged The new formulation consists of evaluating the heat exchanged A B between two connected nodes (a, b) ∈ × (corresponding to A B A ∫ A Ψa hct (𝜃 − 𝜃 ) ds in Eq. (3) and ∫ B Ψb hct (𝜃 − 𝜃 B ) ds in Eq. (4)) by a unique quantity Q(a,b) defined as follows: Q(a,b) = hct (𝜃a − 𝜃b ) S(a,b)
(8)
where 𝜃 a and 𝜃 b denote the temperatures of the two nodes a and b belonging respectively to the surfaces A and B ; hct is the heat exchange coefficient; S(a,b) is the nodal integration surface defined by: ⎛ ⎞ S ⎜ 1ICk (a, b) ∫ B Ψb ds ⎟ S(a,b) = A ⎜ ∑ Ψ ds (9) S ⎜ 1ICk (a, c) ∫ B Ψc ds ⎟⎟ ∫ A a ⎝ c ∈ B ⎠ ∑ The term 1ICk (a, b) ∫ B Ψb ds ∕ c∈ B 1ICk (a, c) ∫ B Ψc ds is a partition coefficient different than 1 when the node a is connected to more than one node b. It is such that: A
(11)
can also be estimated by exchanging the two surfaces A and B but its value will be modified. As for the connectivity procedure (see Remark 1), the finite element formulation proposed is natively not commutative but it always leads to a symmetric equations system.
with
b ∈
a∈
In Eq. (9), the surface A is taken as the integration support. S(a,b)
1 (S + S B ) 2 A
SB =
B
S Ψ ds = S S A ∫ A a
received by B . Therefore, the finite element formulation proposed is symmetric.
The size of the surface of heat exchanged is evaluated by considering the mean values of the parts of each surface A and B belonging to the contact. The mean surface is computed as follows:
a∈
b ∈
∑
Remark 2. The balance of heat exchanged is implicitly ensured because Q(a,b) estimates simultaneously the heat lost locally by A and the one
2.3. Surface of heat exchanged
SA =
A
S(a,b) =
In order to reduce the computing time, S(a,b) can be evaluated only one time once the mesh is build.
Remark 1. The connectivity procedure proposed is not necessary commutative because of the rule 1. Indeed, the subsets Ck are not necessary the same if the surfaces SA and SB are exchanged.
S=
∑
∑
B
1ICk : × → {0, 1} { A B 1 if {a, b} ∈ Ck and (a, b) ∈ × , 1ICk (a, b) ≔ 0 otherwise.
∀a ∈ ,
∑ b ∈
∑ B
1ICk (a, b) ∫ B Ψb ds
c ∈
B
1ICk (a, c) ∫ B Ψc ds
=1
(10)
Fig. 5. 2D reference temperature distribution (◦ C). 15
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 8. 3D temperature distribution (◦ C) between 2 meshes of 6 mm thickness made of P1 (left) and P2 (right) finite elements.
Fig. 6. 2D temperature distribution (◦ C) with P1 elements.
new formulation developed in this work and the technique proposed by Feulvarch et al. [18]. The results are satisfactory for the new technique whereas the temperature distribution is unrealistic for the second method (see Fig. 6). The contact seems to be too stiff because the temperature is constant on all the contact surfaces. Different 2D tests have been carried out on Fig. 7 by using higher order P2 and Q2 elements for the mesh of the part located on the right. One can note that the results are very satisfactory despite the high difference in the mesh sizes. Fig. 8 shows the 3D temperature distribution in the case a contact between a mesh made of P1 elements and a second one made of larger P2 elements. This calculation shows once again the relevance of the formulation proposed in this paper for 3D applications.
Fig. 9. Sketch of the thermal problem studied.
3.2. Convergence test This application is similar to the one of the previous application with one added hole as shown in Fig. 9 for having a non-linear spatial evolution of the temperature. We propose to test the convergence of the numerical solution 𝜃 h obtained with P1 finite elements to a reference solution 𝜃 plotted in Fig. 10. The field 𝜃 was built by means of a computation on a 2D mesh made of finite elements of type P1 ◦ (h = 2.10−2 mm) without considering contact conditions. In Fig. 11, the order of convergence 𝛼 is identified by means of the slope of the linear regression on the log-log plot of ‖𝜃 − 𝜃 h ‖L2 (Ω) (≈C hα ). Two first tests have been performed by considering a homogeneous element size on the whole mesh. For the first one denoted T1 , 1 and 2 correspond respectively to A and B whereas for the second test T2 , A and B are exchanged. Then, two other tests have been carried out by considering a fix element size on the left part of the mesh equal to the one used to build the reference solution. The element size h only evolves on the right part of the mesh. As above, 1
Fig. 7. 2D temperature distribution (◦ C) for different types of meshes.
Fig. 10. 2D reference temperature distribution (◦ C). 16
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18
Fig. 11. Convergence plot of the numerical solution 𝜃 h .
and 2 correspond respectively to A and B for the test T3 whereas they are exchanged for the test T4 . One can note that a first order accuracy is obtained for the tests T1 and T2 corresponding to a homogeneous evolution of the element size h. Moreover, the permutation of the contact surfaces has no significant effect on the quality of the solution and the convergence rate. For the tests T3 and T4 , it is difficult to analyse the value of the convergence
Fig. 13. 3D temperature distribution (o C) between 2 meshes of 6 mm thickness made of P1 (left) and P2 (right) finite elements.
rate because there is not a homogeneous evolution of the element size. The interest of these tests is to show that the permutation of the contact surfaces does not affect the quality and the convergence of the solution in the case of non-matching meshes. 3.3. Contact with thermal contact resistance In this section, we consider a contact heat transfer coefficient hct of 103 W∕(m2 .◦ C) in Eq. (8). Then, the temperature field is not continuous from one side to the other side of the contact. The 2D temperature distribution is plotted in Fig. 12 for the contact between a mesh made of P1 elements on the left and a mesh made of P2 or Q2 elements on the right. The corresponding distribution in 3D is shown in Fig. 13. Once again, the results are very satisfactory compared to the reference solution in Fig. 12. 4. Conclusion The aim of this work was to propose a new formulation for the numerical treatment of thermal contact problems without any constraint on the discretization of the contact surfaces and on the order of the meshes in contact. The results are very satisfactory for all the examples and they show the relevance of the formulation developed for thermal contact problems. References [1] P. Wriggers, G. Zavaris, Thermomechanical contact - a rigorous but simple numerical approach, Comput. Struct. 46 (1993) 47–53. [2] L. Johansson, A. Klarbring, Thermoelastic frictional contact problems: modelling, finite element approximation and numerical realization, Comput. Methods Appl. Mech. Eng. 105 (1993) 181–210. [3] V. Robin, A. Sanchez, T. Dupuy, J. Soigneux, J.M. Bergheau, Numerical simulation of spot welding with special attention to contact conditions, Graz-Seggau, Austria, October 2001, in: H. Cerjak (Ed.), Mathematical Modelling of Weld Phenomena 6, The Institute of Materials, London, 2002. [4] D. Pantuso, K.-J. Bathe, P.A. Bouzinov, A finite element procedure for the analysis of thermo-mechanical solids in contact, Comput. Struct. 75 (6) (2000) 551–573. [5] R.N. Raoelison, A. Fuentes, C. Pouvreau, Ph Rogeon, P. Carre, F. Dechalotte, Modeling and numerical simulation of the resistance spot welding of zinc coated steel sheets using rounded tip electrode: analysis of required conditions, Appl. Math. Model. 38 (2014) 2505–2521. [6] Mohsen Loh-Mousavi, Mehdi AhmadiRad, Tomoyoshi Maeno, Denis J. Politis, LiLiang Wang, Coupled thermal-electrical finite element analysis of electrical resistance heating in hot stamping of ultra-high strength steel tubes, Procedia Manuf. 15 (2018) 1047–1054.
Fig. 12. 2D temperature distribution (◦ C) for different types of meshes. 17
E. Feulvarch et al.
Finite Elements in Analysis and Design 162 (2019) 13–18 [16] N. Stromber, Finite element treatment of two-dimensional thermoelastic wear problems, Comput. Methods Appl. Mech. Eng. 177 (3) (1999) 441–455. [17] P. Wriggers, C. Miehe, Contact constraints within coupled thermomechanical analysis - a finite element model, Comput. Methods Appl. Mech. Eng. 113 (1994) 301–319. [18] E. Feulvarch, V. Robin, J.M. Bergheau, Resistance spot welding simulation: a general finite element formulation of electrothermal contact conditions, J. Mater. Process. Technol. 153154 (2004) 436–441. [19] H. Xing, A. Makinouchi, Three dimensional finite element modeling of thermomechanical frictional contact between finite deformation bodies using R-minimum strategy, Comput. Methods Appl. Mech. Eng. 191 (3738) (2002) 4193–4214. [20] ANSYS, Users Manual, 2018. [21] ESI Group, Users Manual, 2018. [22] David M. Trujillo, Chris G. Pappoff, A general thermal contact resistance nite element, Finite Elem. Anal. Des. 38 (2002) 263–276. [23] Seungmin Jin, Dongwoo Sohn, Im Seyoung, Node-to-node scheme for three-dimensional contact mechanics using polyhedral type variable-node elements, Comput. Methods Appl. Mech. Eng. 304 (2016) 217–242. [24] A. Rieger, P. Wriggers, Adaptive methods for thermomechanical coupled contact problems, Int. J. Numer. Methods Eng. 59 (2004) 871–894. [25] E. Feulvarch, J.M. Bergheau, J.B. Leblond, An implicit finite element algorithm for the simulation of diffusion with phase changes in solids, Int. J. Numer. Methods Eng. 78 (2009) 1492–1512.
[7] C. Maniere, A. Pavia, L. Durand, G. Chevallier, K. Afanga, C. Estournes, Finite-element modeling of the electro-thermal contacts in the spark plasma sintering process, J. Eur. Ceram. Soc. 36 (2016) 741–748. [8] C. Agelet de Saracibar, Numerical analysis of coupled thermomechanical frictional contact problems. computational model and applications, Arch. Comput. Methods Eng. 5 (3) (1998) 243–301. [9] T.A. Laursen, On the development of thermodynamically consistent algorithms for thermomechanical frictional contact, Comput. Methods Appl. Mech. Eng. 177 (3) (1999) 273–287. [10] V.G. Oancea, T.A. Laursen, A finite element formulation of thermomechanical rate-dependent frictional sliding, Int. J. Numer. Methods Eng. 40 (23) (1997) 4275–4311. [11] M. Dittmann, M. Franke, I. Temizer, C. Hesch, Isogeometric analysis and thermomechanical mortar contact problems, Comput. Methods Appl. Mech. Eng. 274 (2014) 192–212. [12] S. Hueber, B.I. Wohlmuth, Thermo-mechanical contact problems on non-matching meshes, Comput. Methods Appl. Mech. Eng. 198 (15) (2009) 1338–1350. [13] J. Molinari, M. Ortiz, R. Radovitzky, E. Repetto, Finite-element modeling of dry sliding wear in metals, Eng. Comput. 18 (2001) 592–610. [14] A. Seitz, W.A. Wall, A. Popp, A computational approach for thermo-elasto-plastic frictional contact based on a monolithic formulation using non-smooth nonlinear complementarity functions, Adv. Model. Simul. Eng. Sci 5 (1) (2018) 5. [15] I. Temizer, Sliding friction across the scales: thermomechanical interactions and dissipation partitioning, J. Mech. Phys. Solids 89 (2016) 126–148.
18