A Simple FEM Formulation Applied to Nonlinear Problems of Impact with Thermomechanical Coupling

The thermal effects of problems involving deformable structures are essential to describe the behavior of materials in feasible terms. Verifying the transformation of mechanical energy into heat it is possible to predict the modifications of mechanical properties of materials due to its temperature changes. The current paper presents the numerical development of a finite element method suitable for nonlinear structures coupled with thermomechanical behavior; including impact problems. A simple and effective alternative formulation is presented, called FEM positional, to deal with the dynamic nonlinear systems. The developed numerical is based on the minimum potential energy written in terms of nodal positions instead of displacements. The effects of geometrical, material and thermal nonlinearities are considered. The thermodynamically consistent formulation is based on the laws of thermodynamics and the Helmholtz free-energy, used to describe the thermoelastic and the thermoplastic behaviors. The coupled thermomechanical model can result in secondary effects that cause redistributions of internal efforts, depending on the history of deformation and material properties. The numerical results of the proposed formulation are compared with examples found in the literature.


INTRODUCTION
The current paper presents the study of thermomechanical response for frictionless impact problems, considering the nonlinear behavior through a simple formulation based on FEM, called positional.

Latin American Journal of Solids and Structures 14 (2017) 2439-2462
The formulation is described in Total Lagrangian reference using as unknowns the nodal positions instead of the displacement. This formulation shows a simple language in respect to nonlinear geometric approach where the main advantage is the absence of co-rotational axes (the formulation can be carried on directly with global axes only) Carrazedo and Coda, 2010;Coda and Paccola, 2011;Coda et al., 2013).
Duhamel in 1837 and Neumann in 1885 identified the modification need of the mechanical models to consider the temperature differences to accurately represent the behavior of the materials (Sherief et al., 2004). The thermomechanical theory investigates the interaction between temperature and structural strains. In many cases, the influence of the strains in relation to the thermal field may be neglected. This results in thermomechanical uncoupled system, for which only the effect of temperature modifies the strains field. The experimental studies and theoretical models of thermomechanical have advanced simultaneously. The use of more sophisticated experiments results in more accurate and efficient mathematical models.
The thermoelasticity and thermoplasticity are defined as part of the thermomechanical theory that determines the behavior of elastic and plastic bodies, respectively, submitted to thermal and mechanical efforts. The first studies were devoted to static problems. Danilovskaya (1950) was the first researcher that tried to include the effects of inertia in thermoelastic transient problems. These problems presented simple nature, but the solution explained the process of transmission of the thermal stresses. The theory of uncoupled thermoelasticity may not be satisfactory for some transient problems, because experimental observations that shown the thermal field modified by the strains field. Several experiments emphasize the influence of strains on thermal field (Rittel, 1998;Benzerga et al., 2005;Stanley, 2008). Biot (1956) developed the classical theory of coupled thermoelasticity. The equations of elasticity and heat conduction are interdependent. The equations used to describe this theory are a combination between the elasticity theory and the laws of thermodynamics. The equations that represent the heat conduction are parabolic and heat propagation velocity in an elastic system is infinite. The theory of generalized thermoelasticity came to deal with these inconsistencies observed in experimental results.
The general theory of thermoelasticity has been modified over the years. The first theory proposed by Lord and Shulman (1967) considers that the heat conduction law is altered to consider the heat flux and its rate. It was considered a time parameter relaxation to ensure finite velocity of propagation of heat. The second theory developed by Green and Lindsay (1972) considers that the constitutive equations are modified by insertion of two time parameters relaxation. Green and Naghdi (1993) presented a new theory of generalized thermoelasticity based on the balance of energy and entropy, where energy dissipation is not allowed. Subsequently, Green and Naghdi (1995) presented the most general case for the theory of generalized thermoelasticity, where energy dissipation was taken into account.
The theories of thermoelasticity cannot accurately describe the behavior of elastoplastic bodies. Thus, Dillon Jr. (1963) developed the theory of thermoplasticity with the purpose of represent the evolution of plastic strains. Since then, heat conduction problems with elastoplastic characteristics have been studied extensively and various exact and numerical solutions can be found in the literature, such as McKnight and Sobel (1977), Dargush and Banerjee (1991), Hakansson et al. (2005) and Carrazedo and Coda (2010).
Although the solution of a heat conduction problem does not present major complexity for most materials, the modeling strategy of thermal and mechanical interdependence continues to be a challenge for nonlinear problems. The study nonlinear thermomechanical behavior problems remains subject of many researches in different fields, such as Rajagopal (1995), Rajagopal et al. (1996), Canadija and Brnic (2009), Canadija and Brnic (2010), Ozakin and Yavari (2010) and Yavari and Goriely (2013).

THE COUPLED THERMOMECHANICAL
The thermoelasticity has been used in several areas, especially for the application of static problems (Copetti, 1999;Shahani and Nabavi, 2007) and dynamic problems (Chen and Dargush, 1995;Norris, 2006;Shahani and Bashusqeh, 2014). Other analysis methods have been adopted to deal with coupled thermoelastic problems. As an example, Soler and Brull (1965) used perturbation techniques and more recently Lychev et al. (2010) determined a closed form solution by an expansion of functions generated by the heat conduction and equations of motion.
Several studies described the thermoelastic behavior of orthotropic materials, such as Lu and Pister (1975), Vujosevic and Lubarda (2002) e Lubarda (2004). Another field of research is the study of anisotropic materials (Deb et al., 1991;Li, 1992;Clayton, 2013;Mahmoud et al., 2015). Numerical approximations for thermoelastic equations are commonly found using the FEM. Models for transient thermoelastic FEM are developed and compared with analytical solutions, as presented in Nickell and Sackman (1968) and Ting and Chen (1982). Rand and Givoli (1995) developed a dynamic model thermoelastic for FEM. Formulations of finite elements for thermoelastic damping are obtained from an irreversible entropy flux due to the heat fluxes caused by the variations of volumetric stresses, as presented in Serra and Bonaldi (2009).
The determination of the temperature field in an elastic body is obtained by solving the differential equation of heat conduction, subject to initial and boundary conditions. The thermomechanical behavior of an elastic and heat conducting body is described by the heat conduction equation (Eq. (1)) and local dynamic equilibrium (Eq. (2)), which are the main equations of the theory of coupled thermoelasticity.
, 0 1 2 1 2 n q a qe r q r n where G is the transverse modulus of elasticity, n is the Poisson coefficient, r is the density, e are the deformations, d ij is the Kronecker delta and i F are the external forces applied. The Eqs.
(1) and (2) also show the coefficient of thermal expansion of material ( a ), specific heat ( e c ), coefficient of thermal conductivity ( k ), internal heat source ( R ), reference temperature ( 0 q ) and temperature variation ( q ).

Latin American Journal of Solids and Structures 14 (2017) 2439-2462
The heat variation is defined by the heat flux and the heat generated internally. The equation of thermoplasticity firm in the first and second law of thermodynamics. In the form of the inequality of Clausius-Duhem, one has: where u is the internal energy, s ij are the axial stress, q is the heat flux and S the entropy.
For all energy transformations there is an increase in the entropy. The total energy involved in the process can not vary, and the entropy increase is a way to measure the energy that can be converted into mechanical work. The Helmholtz free-energy measures the energy of a system that can be transformed into work. The Eqs. (5) and (6) define, respectively, the Helmholtz free-energy ( Φ ) and its rate (  Φ ) as a function of temperature (T ).

( )
Combining the concept of Helmholtz free-energy with Eqs. (3) and (4), it is possible to obtain the equation that defines the energy conservation equation in terms of the internal dissipation, defined by: Where, r r q s e -- With L representing the internal dissipation of energy.

Internal Variables
In accordance with experimental results it was observed that part of the plastic work was converted into heat (Farren and Taylor, 1925;Taylor and Quinney, 1934). Dillon Jr. (1963) and Perzyna and Sawzcuk (1973) presented the first attempts to develop constitutive models considering the interaction between the plastic work and thermal effects. Subsequently, formulations have been developed considering with large deformations (Lemonds and Needleman, 1986;Simo and Miehe, 1992;Canadija and Brnic, 2004). Simo and Miehe (1992) present an approach to thermoplasticity analysis, considering a thermodynamically consistent formulation of the problem coupled and detailing the performance and numerical aspects involving the implementation by FEM. Continuing the development of energy conservation equation is necessary to develop the internal dissipation term. Two variables are Latin American Journal of Solids and Structures 14 (2017) 2439-2462 introduced: the plastic deformation ( e p ij ) and the hardening variable ( x i ). The decomposition of the strain tensor is defined additively, where the total strain is the result of the sum of elastic and plastic parts: where e e ij is the elastic deformation. From Eq. (9), the rate of Helmholtz free-energy can be expressed by Eq. (10). For reasons of convenience, from that moment the temperature gradient , Substituting Eq. (10) into (4) one has: According to Carrazedo and Coda (2010), by demanding that (11) holds for all admissible thermomechanical process which means that it must satisfy all independent variations of e , q  and  Θ , Eq. (11) is satisfied if and only if: Considering Eqs. (12), (13) and (14) and substituting Eq. (11) into Eq. (7) gives: Deriving the first term of Eq. (15), one has: In the plastic regime, a large amount of plastic mechanical energy is dissipated as heat. However, the plastic work is not completely transformed into thermal energy. Part of this work is dissipated due to interaction between the interfaces of the microstructures that constitute the material. The absorbed energy due to generation and rearrangement of imperfections in the process of plasticity is defined as stored energy of cold work ( E ), given by: Based on the law of Fourier, specific heat settings and stored energy of cold work, the Eq. (16) is rewritten as follows: where, The Eqs. (19), (20) and (21) represent, respectively, the heat generation due to elastic deformation, the dissipated plastic working and the stored energy of cold work.
The dissipation mechanisms are of several studies. Referring the stored energy of cold work, is possible to highlight the theoretical and the experimental studies of Bever et al. (1973), Oliferuk et al. (1993), Rosakis et al. (2000), Mroz and Oliferuk (2002), Rittel et al. (2012) and Kolupaeva and Semenov (2015). These studies emphasize the complexity of characterizing E , because it is dependent on the accumulated plastic deformations. In the absence of information on the microstructural behavior of materials about the variables influence the process and in what quantity, it is convenient to use a constant factor to represent the energy dissipation.
Therefore, it is assumed that the relationship between the plastic working and stored energy of cold work is defined by a constant factor, here denoted by b (Simo and Miehe, 1992;Zhou et al., 1996;Kapoor and Nemat-Nasser, 1998). Thus, we have the following simplification: Combining Eqs. (22) and (18), has the final expression that defines the heat transfer to thermoplastic problems, given by:

The Heat Conduction Discrete Equation
The energy balance presented at Eq. (23) is solved adopting as a reference the initial configuration of the structure, rewritten as: In that m R defines the heat generated due to mechanical deformations, expressed as: 0 1 2 1 2 n aq e bs e n The heat problem is solved before the mechanical problem. Therefore, employ the heat source from the previous time step t . Thus, for the current time step + D t t, the expression (24) is rewritten as follows: The numerical procedure starts by substituting the function q by a finite element approximation as: where q i and f i defined temperature and shape function in node i , respectively.
To approximate Eq. (26) the method of weighted residues is adopted here. Specifically the Galerkin method, thus: where, where V is the volume, W the work, j w are the weighting functions and f j represent arbitrary constants related to nodes j of the elements.
Through Eqs. (27), (28) and (29), for any value of j w and q i being constant, one has: Manipulating Eq. (30), obtain a similar expression, defined by: The application of the divergence theorem in the third term of Eq. (31) gives: The Eq. (32) expresses the thermodynamic forces applied on the boundary. Therefore, combining Eqs. (31) and (32) one has: , , Developing the volume integrals, Eq. (33) can be written in a matrix form as: where, , ,

POSITIONAL FINITE ELEMENT METHOD
In particular, the study developed uses an alternative formulation, called positional FEM, considers node coordinate positions as variables instead of displacements. The positional formulation is classified as Total Lagrange Formulation (Wong and Tin-Loi, 1990). Although recent, there are many studies that use positional formulation.  present a nonlinear dynamic analysis of one-dimensional structures using Newmark's temporal integration algorithm. Coda and Paccola (2008) study the geometric nonlinear analysis of shells with thickness variation and use of curved elements. Carrazedo and Coda (2010) applied the positional formulation in the study of the thermomechanical coupling in nonlinear problems of impact between trusses and rigid obstacle, through the positional finite element method. Greco et al. (2012) compares the numerical results of the positional and co-rotational formulation for truss problems. It is also worth mentioning the following studies that use the proposed formulation in non-linear problems: , Greco et al. (2013), Reis and Coda (2014), Sampaio et al. (2015) and Siqueira and Coda (2016).
The positional formulation uses the Lagrangian description that describes the kinematics of the deformation in terms of a coordinate system, fixed in space. The principle of minimum potential energy, applying a total Lagrangian description, is applied (Greco et al., 2013). The total potential energy (П) is written by: where e U is the strain energy, C K is the kinetic energy, A K is the loss of energy due to damping and P is the potential energy of concentrated forces applied to the body. The kinetic energy is zero for static problems.
According to Eqs. (39) and (40), the total deformation energy is defined by the integral of the specific strain energy ( e u ) over the initial volume and potential energy of the applied forces is expressed as a function of the external forces applied ( i F ) and position ( i X ). The index i refers to the degree of freedom that are associated forces and positions.
The kinetic energy is given by: Substituting Eqs. (39), (40) and (41) into (38) one has: where, where m c is the damping coefficient. This energy function can be evaluated substituting the exact position field for an approximate non-dimensional field ( x ).
Thus, the position of dynamic equilibrium is defined using the minimum potential energy theorem, by differentiating Eq. (44) regarding the generic nodal position s X , resulting in: Substituting Eq. (43) into (45) one has: The equilibrium Eq. (28) is nonlinear regarding i X . Thus, to dissipate the residual forces is inevitable to use a numerical strategy to solve this problem. In the current work, the Newton-Raphson procedure was used to reach the equilibrium. In order to solve it, a Taylor expansion regarding i X is used as follows: Neglecting higher-order errors 2 O one has: where, , 0 In Eq. (49), the first term represents the Hessian matrix. Thus, the dynamic nonlinear problem is achieved by combining the iterative Newton-Raphson procedure with a temporal integration algorithm.
The total deformation ( e ) of the problem is given by the sum of the elastic ( e e ), plastic ( e p ) and thermal ( q e ) parcels according to the equation: where E is the elastic longitudinal modulus.
In Eq. (51), q e describes the thermal behavior, while the term e p describes the plastic behavior of the element obtained from the constitutive material model.
Here, it will be considered the impact scheme frictionless discussed in Simo et al. (1986), Greco et al. (2004) and Carrazedo and Coda (2010). The so-called scheme of null-penetration condition, having as its basic principle the position limitation of each node of the structure that are impacted. It is used the classic time integration algorithm of Newmark.

Thermal Loading on Rod
In this example, presented by Copetti (2002) is studied the thermal behavior of a one-dimensional bar with thermal loading. The bar was divided into 100 elements of equal size. The author adopts the thermal conductivity, the specific heat and the density of the material equivalent to 1. The temporal discretization is performed with a time interval equal to 0.0001.
The Figure 1 shows the thermal load along the bar, described by Equation: ( ) 10 cos(2 ) The temperature and displacement are constrained at the position 0 = x : The Figures 1 and 2 give a summary of the temperature variations over time. The results obtained are similar to the Copetti (2002).  As shown in Figure 1, because the boundary conditions the temperature of all nodes tend to remain in equilibrium when 4 = t . It should be noted that over time the thermal load is dissipated between the nodes. Figure 2 illustrates the path of the thermal balance for different positions of the bar, where they tend to a common point ( 10 q = ).

Temperature Evolution in a Rod
This example was originally presented by Kamlah and Haupt (1998), and subsequently by Carrazedo and Coda (2010), which consist in a cylindrical rod of 10 cm length under elastoplastic loading and a reference temperature of 293 K , under the following initial and boundary conditions: The following material properties were set: With the aid of the mechanical model response, this process can be divided further in periods of elastic and elastoplastic loading: Confronting the results obtained with the work mentioned above. The different elastic and elastoplastic periods in the history of loading can be easily identified in Figure 3, which illustrates the thermomechanical heat source through time. It is noted that the elastoplastic hardening model has a heat source larger than the perfect elastoplastic model and therefore more relevant temperature changes. The Figure 4 shows the temperature variation in the center of rod.  The Figure 5 shows the temperature variation along the rod for different time instants. Changes in heating and cooling phases are defined by the state of the body changes. As expected, it is noted that the elastoplastic model with strain hardening kinematic presents higher temperature variations than the perfect elastoplastic model, due to the levels of stresses in the kinematic model are higher.

Impact of a Bar on a Rigid Wall
This example, approached initially by Armero and Petocz (1998), consists of uniaxial impact from a bar (with a constant velocity) and rigid wall, as shown in Figure 6. The problem is modified to consider the thermomechanical effects. The geometrical and material characteristics of both elements are given by: , 0.17 a = and 0.05 d = . Only the extremity node of the bar is impacted. It is investigated the behavior of contact forces, velocities, stresses and temperature changes of the bar. In this example it is investigated only geometric nonlinearity. The time discretization is done to 250 time steps of 0.01. The bar is discretized with 20 finite elements of the same dimension.
The Figures 7 and 8 describe respectively the velocity and contact force, comparing the analytical response to the mechanical and thermomechanical numeric responses. In both figures, before the impact, mechanical and thermomechanical response are equal and the structure does not show deformations. The velocity results of mechanical and thermomechanical models of the impactor node are practically identical up to the instant of 1.93. After this time instant, the difference between the models increases with time. As shown in Figure 8, the heat generation resulting in higher contact forces to mechanical response, and consequently the contact time is reduced.
The Figure 9 shows the temperature field to the undeformed configuration of the bar. It is considered the coupled and the uncoupled problem, respectively. It is observed that the thermomechanical coupling causes a secondary effect and temperature variations were more relevant when considering that the temperature changes generate deformations. The bar has maximum temperatures between times instants of 1.0 and 1.5, approximatively. The bar begins to cool quickly after reflection. In the surface of Figure 10 it is possible to compare the stresses on the bar. It is observed a redistribution of stresses caused by thermomechanical coupling. In mechanical problem the stress varies between -0.50 and 0.20. While the problem coupled, has higher stress levels in compression and traction varying from -0.60 to 0.30. Consequently, the thermomechanical coupling increases the contact force and reduces the contact time.

Unidirectional Impact Between Two Bars
In this example it is studied the case of impact between two identical bars with the same initial velocity (Figure 11), however, moving in opposite senses. This example is present in the studies by Carpenter et al. (1991). The geometrical and physical characteristics of both elements are given by:  In Figure 13, verifies the plasticizing effect of heat generation at 9.5 = x in. It is observed that the evolution of temperature changes is very dependent on the accumulation of irreversible deformations. The Figures 14, 15 and 16 compare the temperature changes obtained with different coefficient of thermal expansion. It is noted that smaller coefficients of thermal expansion result in smaller temperature changes and small difference between the uncoupled and coupled responses. The thermomechanical coupled generates a secondary effect which causes larger variations in temperature and changes in the structural behavior.
The results presented in Figure 16 show high levels of temperature. In this case, due to redistribution efforts, the field of temperatures for coupled problem varies between 6.3 º -F and 69.4 ºF , while the uncoupled problem varies between 2.8 º -F and 29.7 ºF . The larger coefficients of thermal expansion result in larger changes in temperature and contact forces and, consequently, the lower the contact time.

Impact Between a Plane Truss and a Rigid Wall
This example is the impact between a circular truss and rigid wall ( Figure 17). The structure has 264 bars and 97 nodes. The time discretization is done through the Newmark method with 0.00001 D = t s. The geometrical and material data of both elements are given by:  Five nodes are impacted ( Figure 17). The Figure 18 shows the evolution of temperature changes of the nodes impacting over time. The temperature of these nodes begin to increase rapidly from the moment of impact. It is observed a small decrease in the variations of temperatures from a certain moment, characterized by Gough-Joule effect and the dissipation of temperatures between nodes. The Figure 19 shows for different time, the temperature distribution of the structure in the deformed configuration.
Therefore, in Figures 18 and 19 emphasize the importance of considering the thermomechanical behavior in impact problems because the high strain rates can cause changes in the structural configuration.

CONCLUSIONS
A simple and effective alternative formulation to deal with the dynamic nonlinear systems, written in relation to nodal positions has been successfully applied to coupled thermomechanical problems. The formulation includes a complete treatment of the analysis of elastoplastic materials for simple configurations.
Through the numerical examples it is possible to observe the difference between the mechanical and the thermomechanical responses. Therefore, it is necessary to consider the interaction between thermal and mechanical behavior, which can cause redistributions of stress. The reference temperature is a determining factor in thermomechanical analysis, and high reference temperatures imply high temperature variations. The study of structures with elastoplastic behavior is important because depending on the history of deformation and accumulation of irreversible deformations, has a significant part in the generation of heat.
It emphasizes the consideration of thermomechanical coupling in impact problems, because the interaction between the mechanical and thermal response cause significant secondary effects due to high strain rates by modifying the structural behavior. Therefore, depending on the material, loading and initial conditions, the thermomechanical coupling can result in predominant contributions in structural response.