Abstract

A new monolithic solution scheme for thermo-elasto-plasticity and thermo-elasto-plastic frictional contact with finite deformations and finite strains is presented. A key feature is the reformulation of all involved inequality constraints, namely those of Hill’s orthotropic yield criterion as well as the normal and tangential contact constraints, in terms of non-smooth nonlinear complementarity functions. Using a consistent linearization, this system of equations can be solved with a non-smooth variant of Newton’s method. A quadrature point-wise decoupled plastic constraint enforcement and the use of so-called dual basis functions in the mortar contact formulation allow for a condensation of all additionally introduced variables, thus resulting in an efficient formulation that contains discrete displacement and temperature degrees of freedom only, while, at the same time, an exact constraint enforcement is assured. Numerical examples from thermo-plasticity, thermo-elastic frictional contact and thermo-elasto-plastic frictional contact demonstrate the wide range of applications covered by the presented method.

Keywords

Introduction

In many engineering applications frictional contact and elasto-plastic material behavior come hand in hand. Just one class of typical well-known examples are metal forming and impact/crash analysis, where, at high strain rates, thermal effects need to be taken into account. The thermo-mechanical coupling appears in several forms: firstly and most obviously, there is heat conduction across the contact interface. Secondly, the dissipation of frictional work leads to an additional heating at the contact interface. Thirdly, also plastic work within the structure is transformed to heat. Vice versa, the current temperature may influence the elastic and especially the plastic material response. All this necessitates robust and efficient solution algorithms for fully coupled thermo-elasto-plastic contact problems, which has been an active research topic over the past 25 years. Most contributions, however, focus either on thermo-plasticity or on thermo-mechanical contact, while resorting to relatively simple standard methods for the remaining problem parts.

Early implementations of thermo-elastic contact based on node-to-segment contact formulations in combination with a penalty constraint enforcement can be found in [1–7]. Within the last decade, more sophisticated variationally consistent contact discretizations based on the mortar method have been developed and applied to thermo-mechanical contact in [8–12]. In addition, those algorithms satisfy the contact constraints exactly (at least in a weak sense) by using either Lagrange multipliers or an augmented Lagrangian functional instead of a simple penalty approach. Due to an easier implementation and other benefits like symmetric operators, most of the cited works above employ some sort of partitioned solution scheme for solving the structural problem (at constant temperature) and thermal problem (at constant displacement) sequentially. In thermo-plasticity, those partitioned schemes based on an isothermal split are only conditionally stable [13]. Only [4, 6, 11, 12] employ monolithic solution schemes, which solve for displacements and temperatures simultaneously. Most developments of advanced computational methods in thermo-mechanical contact are restricted to thermo-elastic effects; coupled thermo-elasto-plastic contact can only be found in [5–7].

Numerical algorithms for finite deformation thermo-plasticity go back to the seminal work by Simo and Miehe [13], which is based on the isothermal radial return mapping algorithm presented in [14, 15]. Both partitioned and monolithic solution approaches are discussed in [13]. Several extensions to this algorithm have been presented later, e.g. a monolithic formulation in principle axes [16] and a variant including temperature-dependent elastic material properties [17]. In a different line of work, a variational formulation of thermo-plasticity has been developed in [18], where the rate of plastic work converted to heat follows from a variational principle instead of being a (constant) material parameter as in [13]. A comparison to experimental results is presented in [19] to support this variational form. We point out that both approaches to determine the plastic dissipation, i.e. [13] and [18], are applicable within the algorithm for thermo-plasticity that will be derived in this manuscript. Besides the mentioned radial return mapping and variational formulations, a different numerical algorithm to isothermal plasticity at finite strains has been developed in [20]. Based on fundamental ideas from [21], the plastic deformation at every quadrature point is introduced as an additional primary variable and the plastic inequality constraints are reformulated as nonlinear complementarity functions. This allows for a constraint violation during the nonlinear solution procedure, i.e. in the pre-asymptotic range of Newton’s method, while ensuring their satisfaction at convergence. As usual in computational plasticity, the material constraints are enforced at each material point independently, such that the additional unknowns can be condensed directly at quadrature point level. It could be shown in [20] that due to this less restrictive formulation, a higher robustness can be achieved, which allows for larger time or load steps.

The present paper now aims at developing a monolithic solution scheme for the thermo-elasto-plastic frictional contact problem based on a new approach. Mortar finite element methods with dual Lagrange multipliers are applied for the contact treatment using nonlinear complementarity functions to deal with both the inequality constraints arising from frictional contact as well as plasticity in a unified manner. This bears novelty both for the numerical formulation of anisotropic thermo-plasticity within the bulk material as well as for the fully nonlinear thermo-mechanical contact formulation at the interface. Furthermore, full compatibility of the algorithms for thermo-plasticity and thermo-mechanical contact is demonstrated. Concerning plasticity, an extension of [20] to coupled thermo-plasticity within a monolithic solution framework is presented. Similar to the isothermal case, the use of Gauss-point-wise decoupled plastic deformation allows for a condensation of the additionally introduced plastic unknowns, where now also thermo-mechanical coupling terms have to be accounted for. The novel thermo-mechanical contact formulation represents a fully nonlinear extension of [12] including a consistent linearization with respect to both the displacement and temperature unknowns. Moreover, the use of dual Lagrange multipliers within a mortar contact formulation enables the trivial condensation of the discrete contact Lagrange multipliers such that the final linearized system to be solved consists of displacement and temperature degrees of freedom only. Our new thermo-mechanical contact formulation is applicable for both classical finite elements based on Lagrange polynomial basis functions as well as isogeometric analysis using NURBS basis functions, for which an appropriate dual basis has recently been proposed in [22]. Owing to the variational basis of the mortar method, the thermo-mechanical contact patch test on non-matching discretizations is satisfied exactly and optimal convergence rates are achieved. Since this paper touches on various topics, and not every reader may be familiar with every single topic, we try to give a self-contained and rather detailed description of the different sub-problems and solution approaches. Even though this requires to some extent the repetition of methods developed elsewhere, we hope to thereby make the article and its novelties amenable to a broader audience.

The remainder of this paper is outlined as follows: “Thermo-plasticity in the bulk continuum” section contains the treatment of thermo-plasticity within the bulk structure from the underlying continuum mechanics to the final discrete system. “Thermo-mechanical contact” section then incorporates thermo-mechanical contact, again starting from a continuum description and closing with the linearized system that needs to be solved in each Newton iteration step. Finally, several challenging numerical examples in “Numerical results” section demonstrate the high accuracy and robustness that can be achieved for benchmark tests and more complex applications in thermo-elasto-plasticity, thermo-elastic contact and thermo-elasto-plastic contact.

Thermo-plasticity in the bulk continuum

Before introducing the new algorithmic treatment of thermo-plasticity in “Solution algorithm using nonlinear complementarity functions” section, we summarize the well-known constitutive relations of thermo-mechanics in “Continuum thermo-mechanics and thermo-plasticity” and “Yield criterion and evolution of internal variables” sections, for which more details can be found in the literature, see e.g. [13, 23, 24]. Only for the simplicity of presentation and not due to any particular restrictions of the developed framework, a few assumptions commonly used in e.g. [13, 17] are adopted here and outlined in “Assumptions on the used free energy” section. In “Weak form of the thermo-mechanical problem” and “Spatial discretization of the thermo-mechanical continuum” sections and, the finite element discretization of the problem (see e.g. [24, 25] and many others) is derived, whereas “Time discretization” section outlines a time discretization based on [26, 27]. Finally, in “Solution algorithm using nonlinear complementarity functions” section, those developments are used to construct a novel nonlinear solution procedure for thermo-plasticity using nonlinear complementarity functions. This can be seen as an extension of the authors’ previous work [20], where increased robustness in the isothermal case as compared to classical return mapping algorithms had been demonstrated in the isothermal case.

Continuum thermo-mechanics and thermo-plasticity

Let the closed set \(\Omega \in \mathfrak {R}^3\) be the reference configuration of a body and \({{\mathbf {x}}}\) the current position of a material point \({{\mathbf {X}}}\in \Omega \) at time t defined by a bijective and orientation preserving mapping \(\varphi _t({{\mathbf {X}}})={{\mathbf {x}}}\). Analogously, we define the current configuration \(\Omega _t = \varphi _t(\Omega )\). The surface \(\partial \Omega \) is divided into the Dirichlet and Neumann boundary \(\Gamma ^D_{m}\) and \(\Gamma ^N_m,\,m\in \{u,T\}\) for the displacements u and the temperature T, respectively. In the time interval of interest \(t\in [0,t_E]\), the following initial boundary value problem (IBVP) must hold:

In the IBVP above, \(\rho _0\) denotes the mass density in reference configuration, \({{\mathbf {F}}}=\mathrm {Grad}\varphi _t\) the deformation gradient, \({{\mathbf {P}}}\) the first Piola–Kirchhoff stress tensor, E the internal energy per unit undeformed volume, R an energy source term per unit undeformed volume, \(\eta \) the entropy per unit undeformed volume, \({{\mathbf {Q}}}\) the material heat flux, \(\hat{{{\mathbf {u}}}}\) the prescribed displacements, \(\hat{{{\mathbf {t}}}}_0\) the prescribed Piola–Kirchhoff tractions, \(\hat{T}\) the prescribed temperatures, and \(\hat{Q}_0\) the prescribed surface heat flux per area in reference configuration. If this heat flux also depends on the temperature at the boundary (as in natural convection boundaries), Eq. (9) becomes a Robin-type boundary condition. Finally, \({{\mathbf {u}}}_0\), \(\dot{{{\mathbf {u}}}}_0\) and \(T_0\) define the initial displacements, velocities and temperature at time \(t=0\), respectively. First, we take a closer look at the last term in (5). From the fact that the absolute temperature T is always positive, one can deduce that

where \({{\mathbf {C}}}={{\mathbf {F}}}^\mathrm {T}{{\mathbf {F}}}\) denotes the right Cauchy–Green tensor. With (13) being assured, the Clausius–Duhem inequality (5) reduces to the Clausius–Planck inequality

Next, we turn our attention to the formulation of elasto-plastic kinematics at finite deformations. As basic concept, we use a multiplicative split of the deformation gradient into an elastic part \({{\mathbf {F}}}^e\) and a plastic part \({{\mathbf {F}}}^p\) as initially proposed in [28]:

Additionally, the entropy is decomposed additively into an elastic and a plastic part, i.e. \(\eta = \eta ^e+\eta ^p\). The plastic part \(\eta ^p\) is associated with the entropy of the plastic configuration (e.g. movement of dislocations) and the elastic part follows from lattice distortion. The interested reader is referred to [13] for a more detailed discussion of this aspect. Moreover, we introduce an additional strain-like scalar internal variable \(\alpha ^i\) associated with isotropic hardening. Kinematic hardening may as well be included via an additional tensor valued internal variable, but for the sake of brevity it is not considered in this work; for an application of the plasticity algorithm presented later with kinematic hardening we refer to [20] and for a derivation of the associated continuum thermodynamics we refer to [18]. The internal energy E is defined as a function of the elastic state only, i.e. \(E=E({{\mathbf {F}}}^e,\eta ^e)\). Introducing the Helmholtz free energy \(\Psi = E-T(\eta -\eta ^p)\) one can reformulate the Clausius–Planck inequality and obtains

As usual in finite strain thermo-plasticity, the free energy \(\Psi \) is assumed to be decomposed additively into an elastic energy contribution, an energy contribution due to work hardening and a thermal energy contribution:

Since \({{\mathbf {F}}}\), \(\dot{{{\mathbf {F}}}}\), T and \(\dot{T}\) may take arbitrary values, we obtain the constitutive relations for the first and second Piola–Kirchhoff stress \({{\mathbf {P}}}\) and \({{\mathbf {S}}}\), respectively

with the Mandel stress tensor \(\varvec{\Sigma } = 2{{\mathbf {C}}}^e \nicefrac {\partial \psi ^e}{\partial {{\mathbf {C}}}^e}\) and the plastic velocity gradient \({{\mathbf {L}}}^p = \dot{{{\mathbf {F}}}^p}{{{\mathbf {F}}}^p}^{-1}\). In case of elastic isotropy considered in the remainder of this paper, the Mandel stress tensor \(\varvec{\Sigma }\) becomes symmetric and the plastic velocity gradient \({{\mathbf {L}}}^p\) in (22) can be replaced by its symmetric part \({{\mathbf {D}}}^p = \mathrm {sym}(\dot{{{\mathbf {F}}}^p}{{{\mathbf {F}}}^p}^{-1})\). Since the first two summands do not depend on the temperature directly, they are referred to as mechanical dissipation and the remaining part as thermal dissipation. Finally, the temperature evolution equation is obtained by inserting (15), (21) and (22) into the energy balance (4). After some algebraic manipulations (see [13]), one obtains

In many computational methods for finite deformation thermo-plasticity (e.g. [13, 17, 29]), a simplified form of plastic heat generation is used based on a dissipation factor \(\beta \), sometimes also referred to as Taylor–Quinney factor. This simplification can also be taken here, i.e. the heat sources due to plasticity [(i.e. \(\mathcal {D}_\text {mech}-T\frac{\partial }{\partial T} \mathcal {D}_\text {mech}\) in (23)] are replaced by a fraction of the total plastic power \(\mathcal {P}_{pl}=\varvec{\Sigma }:{{\mathbf {D}}}^p\). Consequently, (23) becomes

In metal plasticity, the dissipation factor is usually assumed to be in the range of \(\beta \in [0.85,1]\). Within the later presented framework for finite deformation thermo-plasticity, both variants of the energy balance (23) and (25) can be implemented with similar computational effort.

Yield criterion and evolution of internal variables

In the previous section, we have developed the elastic and thermal constitutive relations of thermo-elasto-plasticity. We have, however, not yet decided on a specific plasticity model and flow rule to determine the evolution of the internal variables \({{\mathbf {F}}}^p\) and \(\alpha ^i\), with the only restriction being that the evolution equations must obey the dissipation inequality (22). In rate-independent elasto-plasticity, a yield function defines the set of admissible stress states. We will use an orthotropic yield function originally proposed by Hill [30]

which includes the well-known von Mises criterion as a special case of setting \({\mathbb {H}}\) to the deviatoric projection tensor \({\mathbb {P}}^{\mathrm {dev}}\). For general orthotropic materials with the orthogonal axes \({{\mathbf {n}}}_i,\; i\in \{1,2,3\}\), the orthotropic tensor \({\mathbb {H}}\) is defined via the second order structural tensors \({{\mathbf {N}}}_i = {{\mathbf {n}}}_i\otimes {{\mathbf {n}}}_i\) as

Therein, the tensor products of two symmetric second order tensors are defined as \(({{\mathbf {A}}}\otimes {{\mathbf {B}}})_{ijkl} = A_{ij}B_{kl}\) and \(({{\mathbf {A}}}\odot {{\mathbf {B}}})_{ijkl} = \nicefrac {1}{2}(A_{ik}B_{jl}+A_{jk}B_{il})\). The coefficients are determined by the relations of the normal yield stress \(y_{ii}\) in direction of \({{\mathbf {n}}}_i\), and the shear yield stress \(y_{ij},\,i\ne j\) in the \({{\mathbf {n}}}_i-{{\mathbf {n}}}_j\)-plane, respectively, with a reference yield stress \(y_0\):

Assumptions on the used free energy

We want to briefly comment on some simplifying assumptions posed on the used free energy potential (18). Those simplifications are neither more nor less restrictive on the solution approach presented later than they are for classical radial return mapping algorithms for thermo-plasticity. Hence, they can often be found in the literature in a very similar way, for example in [13, 17]. First, it is assumed that the elastic free energy can be decoupled into the following three summands:

As far as the isothermal response is concerned, this split implies a decoupled volumetric and isochoric elastic response, since \({\mathbb {U}}\) only depends on the elastic change of volume determined by the elastic Jacobian determinant \(J^e = \mathrm {det}{{\mathbf {F}}}^e\) and \({\mathbb {W}}\) only depends on the volume preserving part of the elastic right Cauchy–Green tensor \(\bar{\mathbf {C}}^e= {J^e}^{-\nicefrac {2}{3}}\mathbf {C}^e\). For the modeling of metallic materials, this is a widely used assumption that appears in the exact same way in almost every numerical algorithm for finite strain plasticity, see e.g. [14, 15, 24, 31, 32] and many more. The thermo-mechanical coupling is therefore restricted to \({\mathbb {M}}(J^e,T)\). Following [13], we use

which appears to be the logical extension of a linear small strain thermal expansion model to finite deformations. Therein, \(\alpha _T\) is the linear coefficient of thermal expansion and \(T_0\) a given reference temperature. In summary, the elastic potential (33) accounts for thermal expansion, whereas the elastic material properties such as shear and bulk modulus do not depend on the temperature. A temperature dependent bulk modulus may easily be introduced combining \({\mathbb {U}}\) and \({\mathbb {M}}\) without any changes in the subsequent methods; the isochoric strain energy function \({\mathbb {W}}\), however, is assumed to be temperature independent. To this end, the deviatoric part of the Mandel stress \(\mathrm {dev}\varvec{\Sigma }\) does not depend on the temperature, such that the only dependency on the temperature in the yield function (26) is through \(Y^{pl} = Y^{pl}(\alpha ^i,T)\). The elastic heating term \(\frac{\partial {{\mathbf {P}}}:\dot{{{\mathbf {F}}}}}{\partial T}\) in (23) reduces to

and assume all other potentials in (18) to only depend (piece-wise) linearly on the temperature. As a consequence, the specific heat capacity \(C_v\) defined in (23) and further specified in (24) reduces to a constant material parameter.

Weak form of the thermo-mechanical problem

To set the scene for the subsequent finite element discretization, we introduce the weak form of the thermo-mechanical problem. Therefore, appropriate solution and testing spaces \(\mathcal {U}\) and \(\mathcal {V}\) for the displacement field \({{\mathbf {u}}}\) and temperature field T are defined:

The weak form of the coupled thermo-mechanical problem then consists of the balance of linear momentum (2) and heat conduction (23): Find \({{\mathbf {u}}}\in \varvec{\mathcal {U}}_u\) and \(T\in \mathcal {U}_T\), such that:

Additionally, the thermo-plastic constraints (29), (30) and (31) have to be satisfied locally, such that the thermo-mechanical coupling enters the structural equilibrium via the second Piola–Kirchhoff-stress \(\mathbf {S} = \mathbf {S}({{\mathbf {u}}},T,{{\mathbf {F}}}^p,\alpha ^i)\). Vice versa, the thermal constitutive relation in (14) as well as the source terms \(\mathcal {D}_\text {mech}\) in (22) and \(\mathcal {H}^{ep}\) in (23) introduce the coupling of temperature and displacement field as well as the plastic deformation in (42).

Spatial discretization of the thermo-mechanical continuum

The position, displacement and temperature field (as well as their variations) are approximated in space using discrete nodal values (or control point values in case of isogeometric analysis (IGA)) \({{\mathbf {\mathsf{{X}}}}}_j\), \({{\mathbf {\mathsf{{d}}}}}_j \) and \({\mathsf {T}}_j\) and ansatz functions \(N_j\), viz.

The ansatz functions \(N_j\) of node j may be either Lagrange polynomials for classical finite elements or NURBS in the case of isogeometric analysis. The vectors \({{\mathbf {\mathsf{{d}}}}}\) and \({{\mathbf {\mathsf{{T}}}}}\) contain all displacement and temperature degrees of freedom in the approximation, respectively. As usual in finite element methods for plasticity, the plastic constraints (29), (30) and (31) are enforced locally at the quadrature points. The internal variables \({{\mathbf {F}}}^p\) and \(\alpha ^i\) are therefore assumed to be discontinuous and independent at every quadrature point q, denoted as \({{\mathbf {\mathsf{{F}}}}}^p_q\) and \(\upalpha ^i_q\) in the following. Again, the vectors \({{\mathbf {\mathsf{{F}}}}}^p\) and \(\varvec{\upalpha }^i\) represent the union of all discrete values at the quadrature points.

Remark 1

All algorithms presented later are directly applicable to both finite elements and isogeometric analysis, such that no further distinction will be made in the following. For the sake of brevity, no introduction to isogeometric analysis will be given here, since there has been an overwhelming amount of publications on IGA in the past decade, including the monograph [25]. For details on the application of the dual mortar method to isothermal isogeometric contact mechanics, we refer to our recent work [22].

The spatial discretization (43) can now be inserted into the weak form of the balance of linear momentum in (41), while still neglecting the contact contribution for now. The discrete algebraic force equilibrium becomes:

where \({{\mathbf {M}}}_u\) is the (constant) mass matrix. Still, the constraints of elasto-plasticity (29), (30) and (31) have to be applied. The way they are introduced to the discrete system is discussed in more detail in “Solution algorithm using nonlinear complementarity functions” section. It is also due to plasticity (and thermal expansion) that a dependency on the discrete temperatures \({{\mathbf {\mathsf{{T}}}}}\) and the internal variables \({{\mathbf {\mathsf{{F}}}}}^p\) and \(\varvec{\upalpha }^i\) is introduced into the internal force vector via the second Piola–Kirchhoff stress \(\mathbf {S}_q = \mathbf {S}_q({{\mathbf {C}}}_q, \mathrm{T}_q,{{\mathbf {\mathsf{{F}}}}}^p_q,\upalpha ^i_q)\).

The same discretization can be applied to the weak form of the heat conduction equation, which gives

Similar to the structural equilibrium above, the mass matrix \({{\mathbf {M}}}_T\) determining the heat capacity is constant and the internal load vector \({{\mathbf {f}}}_T^\text {int}\) depends linearly on the temperature and nonlinearly on the displacement via Fourier’s law of heat conduction in the finite deformation realm, see (14). The discrete mechanical dissipation vector \({{\mathbf {f}}}^\text {diss}_T\) depends nonlinearly on the displacement and the plastic deformation at every quadrature point according to (22) and (23). Finally, for both the structural and the thermal problem, the external load vectors \({{\mathbf {f}}}_{\{u,T\}}^\text {ext}\) are assumed to be independent of the displacement and temperature field for the sake of simplicity.

Time discretization

To discretize the semi-discrete equilibrium (45) and (47) in time, we apply generalized-\(\alpha \) schemes, which are of second-order accuracy, and can be formulated with the spectral radius in the high frequency limit \(\rho _\infty \) as sole parameter. For structural problems, this method has been presented in [26]. The approximation of discrete velocities \({{\mathbf {\mathsf{{v}}}}}\) and accelerations \({{\mathbf {\mathsf{{a}}}}}\) is based on the Newmark-scheme, viz.

where \(\Delta t\) denotes the time step size of the interval \([{}^{n}{}{t},{}^{n+1}{}{t}]\). The left superscript signifies the approximation at the discrete time \({}^{n}{}{t}\) and \({}^{n+1}{}{t}\), respectively. The discrete equilibrium (45) is then evaluated at a generalized mid-point by introducing the parameters \(\alpha _{u,f}\) and \(\alpha _{u,m}\):

The discrete forces (and accelerations) at the mid-points are eventually interpolated by the forces (and accelerations) at the end of each time step, e.g. \({}^{n+1-\alpha _{u,f}}{}{{{\mathbf {f}}}}^\text {int}_u = (1-\alpha _{u,f}) {}^{n+1}{}{{{\mathbf {f}}}}^\text {int}_u +\alpha _{u,f} {}^{n}{}{{{\mathbf {f}}}}^\text {int}_u\). In [26], an optimal set of parameters is derived in terms of the spectral radius in the high frequency limit \(\rho _{u,\infty }\) as

The generalized-\(\alpha \) method has been extended to systems of first order in time in [27], which will be used for the temporal discretization of the thermal evolution (47). Similar to (48), the temperature rate is approximated by

where the values at the mid-points are again obtained by an appropriate linear combination of the end-point values. An optimal choice of the parameters has been derived in [27] in terms of the spectral radius in the high frequency limit \(\rho _{T,\infty }\) as

Finally, we need a discrete time integration of the evolution equations for the internal plastic variables in (29) and (30). Here, we follow standard techniques in finite strain plasticity, see e.g. [24], namely an exponential map time integration for the plastic deformation gradient and a backward-Euler scheme for the other internal variables, viz.

where the plastic multiplier increment \(\Delta \gamma _q\) must be determined such that the Karush–Kuhn–Tucker conditions (31) are fulfilled at time \({}^{n+1}{}{t}\). The advantage of the exponential map here is that it preserves the plastic incompressibility in the time-discrete setting. This means that if the plastic constitutive equations are such that plastic deformation does not result in a change of volume (i.e. \(\mathrm {det}[{{{\mathbf {F}}}^p}]\equiv 1\;\forall t\)), the argument of the exponential function is traceless, and hence, we also get \(\mathrm {det}[{}^{n}{}{{{\mathbf {\mathsf{{F}}}}}}^p] \equiv 1 \;\forall n\) in the time-discrete setting.

Remark 2

While being relatively easy to implement and fairly robust, the presented generalized-\(\alpha \) time integration schemes are not algorithmically energy conserving. As an alternative, so-called structure preserving time integration schemes based on the (generalized) energy momentum method have been proposed in [33–36] for isothermal nonlinear elasticity. Later, those methods have been extended to isothermal contact [37], elasto-plasticity [38], thermo-elasticity [39–41] and thermo-elastic contact [11]. The combination of the cited works to a structure preserving time integration for a fully coupled thermo-elasto-plastic contact problem is beyond the scope of this paper, but might be a worthwhile topic of future research.

Solution algorithm using nonlinear complementarity functions

A similar approach to the one described in the following has been developed for small strain von Mises plasticity in [21] and extended to finite deformation anisotropic Hill-type plasticity in [20]. We therefore introduce the incremental plastic flow \(\Delta {{\mathbf {\mathsf{{D}}}}}^p_q ={\Delta \gamma _q}\frac{{\mathbb {H}}:{}^{n+1}{}{\varvec{\Sigma }_q}}{\Vert {}^{n+1}{}{\varvec{\Sigma }_q}\Vert _{\mathbb {H}}} \) at each quadrature point q as an additional primary variable and can express the evolution equations (55) solely in terms of this incremental plastic deformation, viz.

at each quadrature point q. Here, the temperature only enters via the temperature dependent yield stress \(Y^{pl}_q\); the (deviatoric part of the) Mandel stress is temperature independent due to the assumed restrictions on the used free energy (33). Temperature dependent elastic material properties (and hence a temperature dependency of \(\mathrm {dev}\varvec{\Sigma }\)) may be included, but would further complicate the derivative \(\nicefrac {\partial C^{pl}_q}{\partial {{\mathbf {\mathsf{{T}}}}}}\) needed later and are therefore omitted here for the time being for the sake of simplicity.

The set of nonlinear equations that needs to be solved in every time step of a thermo-elasto-plastic problem without contact consists of the balance of linear momentum (49) and the heat conduction (53) complemented with the NCP function (59) at every quadrature point. Since this system is semi-smooth, variants of Newton’s method can be applied resulting in the linearized system

where the set \({\mathcal {G}}\) contains all potentially plastifying quadrature points. For brevity, the detailed structure of the involved linearizations are omitted here, since they have either already been presented for the isothermal case in [20] or follow from similar calculations as given therein. The system (60)–(62), however, is of significantly increased size compared to the original number of displacement and temperature degrees of freedom, since at every quadrature point, assuming a symmetric and traceless incremental plastic flow \(\Delta {{\mathbf {\mathsf{{D}}}}}^{p}_q\), additional five unknowns are introduced. The key to eliminating the plastic increment \(\Delta {{\mathbf {\mathsf{{D}}}}}^{pl}_q\) from (60)–(62) is the fact that (62) only contains the discrete plastic increment of one quadrature point as well as the displacement and temperature degrees of freedom belonging to the element containing this quadrature point. Hence, the linearized plastic NCP function (62) can be solved directly for the increments \(\Delta (\Delta {{\mathbf {\mathsf{{D}}}}}^{pl}_q)\) at the level of quadrature points, yielding

This condensation can in turn be inserted into (60) and (61) and results in modified linearizations with respect to the displacements and temperatures in those equations, which will be indicated by a tilde \(\tilde{(\cdot )}\). Specifically, one obtains

After having solved for the displacement and temperature increments \(\Delta {{\mathbf {\mathsf{{d}}}}}\) and \(\Delta {{\mathbf {\mathsf{{T}}}}} \) in the current Newton iteration, the condensation equation in (63) can be used to recover the increments of the plastic deformation \(\Delta (\Delta {{\mathbf {\mathsf{{D}}}}}^{p}_q)\) at each quadrature point q. The matrices containing the linearizations with respect to the plastic deformation hence do not have to be assembled at a global level but the condensation (64)–(66) can be performed locally at every quadrature point. The numerical effort is therefore similar to classical radial return mapping algorithms, for which, when applied to general hyperelastic material, a local system of nonlinear equations needs to be solved at every quadrature point, see e.g. [24]. The involved terms are very similar and basically consist of material models (i.e. derivatives of elastic energies) and kinematic evolution equations of which especially (56) is computationally expensive, since it involves matrix exponentials and (for the linearization) their derivatives. For certain hyperelastic materials, the return mapping schemes can be reduced to solving a scalar nonlinear equation only. Similar modifications might be possible for the presented method, but have not yet been explored. In summary, the treatment of (thermo-) plasticity based on NCP functions does not result in higher computational costs compared to classical return mapping algorithms, while allowing for larger load steps as demonstrated in [20].

Thermo-mechanical contact

Having dealt with the thermo-mechanical coupling within the bulk structure, we turn our focus to thermo-mechanical contact problems. First, the continuum mechanical description is recalled in “Problem setup” section, which, in more detail, can also be found e.g. in [42, 43]. Next, a weak form extending the small strain formulation presented in [12] is derived in “Weak form of the thermo-mechanical contact problem” section. Finally, we present the new mortar finite element formulation for finite deformation thermo-mechanical contact problems in “Mortar finite element discretization of thermo-mechanical contact” section.

Problem setup

We consider a two body finite deformation thermo-mechanical contact problem in three spatial dimensions. Let the closed sets \(\Omega ^{(i)} \subset \mathfrak {R}^3, \, i=1,2\) be the reference configuration of the two bodies. Both bodies are governed by the IBVP described in “Thermo-plasticity in the bulk continuum”, enhanced with the constraints of frictional contact at the potential contact boundary \(\Gamma ^{C(i)}\). The boundary \(\partial \Omega ^{(i)}\) is hence decomposed into three subsets such that \(\partial \Omega ^{(i)} = \Gamma ^{D(i)}_u \cup \Gamma ^{N(i)}_u \cup \Gamma ^{C(i)} = \Gamma ^{D(i)}_T \cup \Gamma ^{N(i)}_T \cup \Gamma ^{C(i)}\) and \(\emptyset = \Gamma ^{D(i)}_u \cap \Gamma ^{N(i)}_u = \Gamma ^{D(i)}_u \cap \Gamma ^{C(i)}=\Gamma ^{N(i)}_u \cap \Gamma ^{C(i)} = \Gamma ^{D(i)}_T \cap \Gamma ^{N(i)}_T = \Gamma ^{D(i)}_T \cap \Gamma ^{C(i)}=\Gamma ^{N(i)}_T \cap \Gamma ^{C(i)}\). Since the focus is on finite deformations, the geometrical contact constraints such as the non-penetration condition have to be satisfied in the current configuration, i.e. they have to be enforced between the current potential contact surfaces \(\gamma ^{C(i)} = \varphi _t(\Gamma ^{C(i)})\). Applying standard nomenclature in contact mechanics, we will refer to \(\gamma ^{C(1)}\) as the slave surface, and to \(\gamma ^{C(2)}\) as the master surface. For a point \({{\mathbf {x}}}^{(1)}\) on the slave surface (in spatial configuration), one can find an associated point \(\hat{{{\mathbf {x}}}}^{(2)}\) on the master surface by projecting \({{\mathbf {x}}}^{(1)}\) along its current outward normal vector \({{\mathbf {n}}}\). With those points at hand, the normal gap \(g_n\) and the relative tangential velocity are defined by

where \(\dot{{{\mathbf {x}}}}^{(1)}\) and \(\dot{\hat{{{\mathbf {x}}}}}^{(2)}\) denote the material velocities of \({{{\mathbf {x}}}}^{(1)}\) and \({\hat{{{\mathbf {x}}}}}^{(2)}\), respectively. The contact traction \({{\mathbf {t}}}_c\) at the interface is decomposed in the same way to obtain the normal contact pressure \(p_n\) and the tangential contact traction \({{\mathbf {t}}}_\tau \), viz.

The mechanical contact constraints can then be formulated in normal direction via the Hertz–Signorini–Moreau condition and in tangential direction by Coulomb’s law of friction on the slave contact surface:

The master side contact traction follows directly from the balance of linear momentum \({{\mathbf {t}}}_c^{(2)} = -{{\mathbf {t}}}_c^{(1)}\). Since we are interested in a fully coupled thermo-mechanical contact problem, we allow for a temperature dependent friction coefficient \(\mu (\theta _c)\) depending on the maximal temperature of the contact surfaces \(\theta _c = \mathrm {max}[T^{(1)}({{\mathbf {x}}}^{(1)}),T^{(2)}(\hat{{{\mathbf {x}}}}^{(2)})]\). Following [42], an exemplary temperature dependency of the friction coefficient is introduced as

via a reference coefficient of friction \(\mu _0\) at the reference temperature \(T_0\) and a damage temperature \(T_d>T_0\). The apparent coefficient of friction decreases monotonically from \(\mu _0\) at \(T_0\) to zero at \(T_d\). The damage temperature \(T_d\) is usually chosen to be the lower melting temperature of the two contacting materials, since at this point, friction is no longer dominated by solid shearing but rather by viscous effects in a thin film of molten material.

Next, the local energy balance at the contact interface is investigated. For simplicity, we assume that the interface has no heat capacity by itself, i.e. the specific internal energy \(e_c\) is zero at all times. Written in rate form one obtains

where \(q^{(i)}_c = {{\mathbf {q}}}^{(i)}{{\mathbf {n}}}^{(i)}\) denotes the spatial (Cauchy) heat flux at the two contact surfaces. A possible choice of the heat fluxes at the contact interface is given in [2, 42] as

Here, the heat transfer parameter \(\gamma ^{(i)}\) has been chosen as a simple linear function of the normal contact pressure. Yet, any nonlinear relation may be employed as long as a zero contact pressure (i.e. no contact) relates to a zero heat flux over the contact interface. Using (74) and (75) and eliminating the contact surface temperature \(T_c\) from the system, the spatial heat fluxes can be stated as

with the temperature jump over the contact interface \([T]=(T^{(1)}-T^{(2)})\) as well as the coefficients \( \beta _c =\frac{\bar{\gamma }^{(1)}\bar{\gamma }^{(2)}}{\bar{\gamma }^{(1)}+\bar{\gamma }^{(2)}} \) and \(\delta _c=\frac{\bar{\gamma }^{(1)}}{\bar{\gamma }^{(1)}+\bar{\gamma }^{(2)}} \) defined according to [2, 12].

Remark 3

The presented thermal conditions at the contact interface are actually thermodynamically consistent, meaning that they obey the first and second law of thermodynamics. To prove this, the conditions above can be reformulated in terms of a dissipation potential at the interface and a subsequent analysis similar to the one given above for the bulk continuum can be conducted, see e.g. [2, 42] for a detailed derivation. Since this derivation would lead to no further insight with regard to the following numerical algorithms, it is skipped here and the reader is referred to the respective literature instead.

Remark 4

The parameters \(\beta _c\) and \(\delta _c\) also have a direct physical interpretation: The product \(\beta _cp_n\) in (76) and (77) determines the heat flux across the contact interface due to the temperature jump, i.e. it describes the contact heat conductance. In the presented description, which follows [2, 42], the heat flux is a linear function of the contact pressure. However, more involved nonlinear relations can be found in the literature. For instance, [43] distinguishes three sources of heat conduction across rough surfaces: conduction through contacting asperities, heat conduction in enclosed gas and radiation. In the formulation presented later, nonlinear models for heat conduction could be employed simply by replacing the product \(\beta _cp_n\) with a nonlinear relation \(\bar{\beta }_c(p_n)\) or, when formulated using the later defined Lagrange multipliers (representing the negative contact traction), by replacing \(\beta _c\lambda ^c_n\) with \(-\bar{\beta }_c(-\lambda ^c_n)\) in (95).

The parameter \(\delta _c\) determines how frictional heat is distributed to the two bodies. In the limit cases \(\delta _c=0\) or \(\delta _c=1\) the entire frictional heat is added to the master or slave side, respectively.

Weak form of the thermo-mechanical contact problem

To prepare the subsequent mortar finite element discretization of the thermo-mechanical contact problem, the weak form of the thermo-structure-interaction problem in (41) and (42) is extended to account for the effects of frictional contact. Therefore, two Lagrange multiplier fields are introduced at the contact interface; a vector-valued Lagrange multiplier to enforce the mechanical contact constraints (70) and (71), which can be identified as the negative slave side contact traction \(\varvec{\lambda }^c = - {{\mathbf {t}}}_c^{(1)}\), and a scalar Lagrange multiplier to enforce the heat flux constraints over the contact interface in (76) and (77), which will be chosen as the slave side heat flux \(\lambda ^\theta = q^{(1)}_c\). The contact Lagrange multiplier is taken from the convex set

Therein, all integrals over the contact surfaces of the two bodies have been transformed into pure slave side integrals by using a suitable mapping \(P_t:\gamma ^{C(1)}\mapsto \gamma ^{C(2)}\) at time t.

Mortar finite element discretization of thermo-mechanical contact

The discrete interpolation of displacements and temperatures at the contact interface follows directly from (43) by restricting the ansatz functions to the element boundary. In the present notation, we will not explicitly distinguish between the ansatz functions in the continuum and at the boundary. Instead, the context will provide the necessary information, i.e. if the quantities are integrated over a volume, the ansatz functions refer to (43) and if integration is performed on the contact surface only the trace space restriction of the ansatz functions are evaluated.

where \(\delta _{ij}\) denotes the Kronecker symbol, i.e. \(\delta _{ij}=1\) if \(i=j\) and \(\delta _{ij}=0\) otherwise. In practice, the easiest way to define those dual basis functions is via an element-wise linear combination of the standard ansatz functions \(N_i\), see e.g. [44], which at least ensures a partition of unity property. For details on their construction, linearization and application to contact mechanics with small and large deformations including friction, the reader is referred, for instance, to [45–50], and for a first application to thermo-elastic contact to [12]. With the discretization of the displacement and temperature field (43) and the discrete Lagrange multipliers (83), we can now approach the contact related parts in (79)–(82) and discretize them adequately in space and time.

First, we consider the contact contribution \(G_u^c\) in the weak balance of linear momentum in (79). Inserting the discretization yields

where the set \(\mathcal {M}\) contains all nodes of the master contact surface and \(P^h_t\) denotes a discrete version of the projection \(P_t\) introduced in (79). The first integral in (85) constitutes a square coupling matrix \({{\mathbf {\mathsf{{D}}}}}_u\), which obviously becomes diagonal when using the dual basis from (84). The second integral constitutes the rectangular coupling matrix \({{\mathbf {\mathsf{{M}}}}}_u\). With those matrix blocks and by re-arranging the vector of displacement degrees of freedom \({{\mathbf {\mathsf{{d}}}}}=[{{\mathbf {\mathsf{{d}}}}}_\mathcal {N},{{\mathbf {\mathsf{{d}}}}}_\mathcal {M},{{\mathbf {\mathsf{{d}}}}}_\mathcal {S}]^\mathrm {T}\), where \({{\mathbf {\mathsf{{d}}}}}_\mathcal {S}\) and \({{\mathbf {\mathsf{{d}}}}}_\mathcal {M}\) contains all degrees of freedom on the slave and master contact surface, respectively, and \({{\mathbf {\mathsf{{d}}}}}_\mathcal {N}\) all other nodal displacements, the discrete contact force is obtained as

The displacement dependency of the coupling matrices is a consequence of the projection \(P^h_t\) and the integration, which has to be performed on the current configuration of the contact surface \(\gamma ^{C(1)h}\). Next, this contact force has to be incorporated in the time-discrete balance of linear momentum (49). This is done in a fully implicit way, so that the space and time discrete equilibrium (49) is complemented with a contact contribution

Remark 5

By choosing a fully implicit time discretization of the contact forces instead of some linear combination of \({}^{n+1}{}{{{\mathbf {f}}}}^c_u\) and \({}^{n}{}{{{\mathbf {f}}}}^c_u\), we guarantee that the contact work in one time step \(W^c = ({}^{n+1}{}{{{\mathbf {\mathsf{{d}}}}}}-{}^{n}{}{{{\mathbf {\mathsf{{d}}}}}})^\mathrm {T}{}^{n+1}{}{{{\mathbf {f}}}}^c_u\) becomes negative (i.e. dissipative) for nodes that come into contact and zero for nodes leaving the active contact set. Following an idea in [51] the dissipated energy can be re-introduced into the system via a velocity update procedure. If, however, the contact force were discretized by a linear combination of two time steps, nodes leaving the active contact set would introduce energy to the discrete system, which one might not be able to compensate via the velocity update procedure.

The next contact contribution to be discretized is the contribution to the heat conduction equation (80). A strict application of the discretization (43) and (83) would give

The first and second part therein stem from the heat conduction over the contact interface and result in the similar coupling matrices \({{\mathbf {\mathsf{{D}}}}}_T\) and \({{\mathbf {\mathsf{{M}}}}}_T\) already used in the structural coupling above (the only difference being the size of the matrices: one entry per node for the thermal part and three entries per node—one for each displacement degree of freedom—in the structural coupling). The last integral is the result of frictional dissipation at the contact interface and, in the stated form gives rise to some complications. First, the employed discrete tangential velocity is not frame indifferent, see e.g. [52, 53]. Moreover, it involves a triple integral at the contact interface, which poses high demands on the quadrature accuracy at the contact interface, especially when going to higher order approximations using Lagrange polynomials or NURBS. Following the work of [12], we seek to reduce the computational cost by appropriate lumping techniques in the discrete representation of (80) and (82). Therefore, we introduce a frame indifferent weighted nodal tangential velocity \(\tilde{{{\mathbf {\mathsf{{v}}}}}}_{\tau j}\) derived from a (discrete) time derivative of the mortar matrices \({{\mathbf {\mathsf{{D}}}}}_u\) and \({{\mathbf {\mathsf{{M}}}}}_u\) as

From a physical point of view, this means that we do not interpolate the velocities and the contact Lagrange multiplier but only the scalar product \(\frac{\tilde{{{\mathbf {\mathsf{{v}}}}}}_{\tau j}\cdot \varvec{\uplambda }^c_{\tau j}}{\int _{\gamma ^{C(1)h}}\Phi _j\,\mathrm {d}\gamma } = \mathcal {D}_{\text {mech},j}^c\), which represents the frictional dissipation power at the slave node. Numerically, (90) implies a lumping of the triple integrals in (88), thus resulting in an integral over the product of two ansatz functions only. The contact contribution to the discrete thermal equilibrium (53) is again interpolated at the mid-point \({}^{n+\alpha _{T,f}}{}{{{\mathbf {f}}}}_T^c = \alpha _{T,f}{}^{n+1}{}{{{\mathbf {f}}}}_T^c + (1-\alpha _{T,f}){}^{n}{}{{{\mathbf {f}}}}_T^c\), such that

where the vector \( \varvec{\mathcal {D}}_{\text {mech}}^c\) contains all entries of the contact dissipation power at all nodes \(\mathcal {D}_{\text {mech},j}^c\).

Still, the interface constraints \(H_u\) and \(H_T\) in (81) and (82) need to be discretized. Owing to the biorthogonality of the dual basis functions for the Lagrange multiplier (84), it is well-known that the discretization of the variational inequality (81) yields decoupled constraints for the nodes [45, 48]. The discrete normal contact constraint becomes

In the fully coupled thermo-mechanical contact problem, the coefficient of friction depends on the temperatures at the two sides of the contact interface, e.g. via (73). In the discrete setting we determine the maximal contact interface temperature \(\theta _{cj} = \mathrm {max}[T^{(1)}_j,T^{(2)}({{\mathbf {\mathsf{{x}}}}}_j^{(1)}\circ {P_{t}^h}^{-1})]\). The algorithmic treatment of those inequality constraints will be the topic of the following “Solution algorithm using nonlinear complementarity functions” section.

Finally, a discrete representation of the thermal interface condition (82) is required. At this point, we again deviate from a strict application of the discretization, as we already did in (90), and end up with the following discrete representation of (82):

Thanks to the lumping procedure, especially for the first integral, the thermal interface condition decouples for the contact nodes j, i.e. we can set \(\uplambda ^\theta \) to zero for all inactive contact nodes. For the contact interface constraints, this decoupling can be achieved in a consistent manner using dual shape functions, such that only local, decoupled constraints have to be solved instead of an inequality constraint coupling all interface nodes. To keep up this advantage in the coupled thermo-mechanical contact problem, the presented lumping is required, see [12] for a more detailed discussion. We want to emphasize that for standard thermo-mechanical mortar methods, e.g. [11], a similar simplification is made implicitly by using a node-wise decoupled active set strategy, although strictly speaking the variational inequality does not not allow for such a decoupled treatment in that case [54, 55].

Solution algorithm using nonlinear complementarity functions

The discrete system derived in the previous section still includes discrete inequality constraints for normal contact (93) and Coulomb friction (94) at all slave nodes \({\mathcal {S}}\). As commonly used in isothermal dual mortar methods, the contact constraints are reformulated as nonlinear complementarity (NCP) functions, which then pose non-smooth equality constraints to the system, see e.g. [45, 47, 54]. A first application to small-strain thermo-elastic contact problems was presented in [12]. However, a fixed-point solution approach for the thermal interface conditions was used, thus sacrificing a quadratic rate of convergence. In this work, the concept of NCP functions will only be briefly reviewed by highlighting the additional complexity resulting from thermo-mechanical coupling. At each slave node, we define the trial values

The fully coupled nonlinear system to be solved for each time step comprises the structural and thermal equilibrium (87) and (92), respectively, the plastic NCP function (59), the contact NCP functions (98) and (99), and, finally, thermal contact interface condition (95). All in all, we obtain

This coupled system of equations is now solved monolithically using (a non-smooth version of) Newton’s method. Therefore, the system is linearized with respect to the discrete displacements, temperatures, Lagrange multipliers and plastic flow variables, which results in

(106)

(107)

(108)

(109)

(110)

(111)

Again, we do not elaborate on the linearizations of the involved contact terms, since many of those have already been presented elsewhere; see e.g. [48, 49] for the linearization of the contact forces and the normal contact constraints and [50] for the isothermal frictional contact. The remaining linearizations, especially those with respect to the discrete temperatures, follow in a rather straightforward manner.

As already presented in “Solution algorithm using nonlinear complementarity functions” section, the discrete plastic flow at each quadrature point can be condensed at a Gauss point level, which also applies to the extended system (106)–(111). Finally, we also seek to eliminate the Lagrange multipliers that were introduced at the contact interface to enforce the normal and frictional constraints as well as the heat production and conduction. Being well-established for isothermal contact in the meantime, the general procedure will only be outlined briefly here. As discussed above, the displacement and temperature degrees of freedom will be split into different sets, where \(\mathcal {S}\) and \(\mathcal {M}\) contain the degrees of freedom on the slave and master side, respectively, and \(\mathcal {N}\) all other nodes. Then, we can split the rows of the linearized system (106)–(111) into those sets and obtain a system that can be written in the following matrix form:

By splitting the slave set \(\mathcal {S}\) further into an active part \(\mathcal {A}\), with nodes j where \(\uplambda _{nj}^c-c_n\tilde{g}_{nj}>0\), and an inactive part \(\mathcal {I}\), we can omit the constraints for the inactive nodes, since they only yield \(\varvec{\uplambda }^{\{c,\theta \}} = 0\) for those nodes. The effective stiffness blocks \(\tilde{{{\mathbf {K}}}}\) therein contain the mass matrices \({{\mathbf {M}}}_{\{u,T\}}\) as well as the condensed plastic constraints (64)–(66). The remaining blocks \({{\mathbf {S}}}\)–\({{\mathbf {Z}}}\) can be identified by comparing (112) with the linearized system (106)–(111) and are not specified in detail here. What is crucial for the condensation of the Lagrange multipliers is the fact that the coupling matrices \({{\mathbf {\mathsf{{D}}}}}_{\{u,T\}}\) are of square and diagonal shape due to the biorthogonality property of the dual shape functions. Hence, the Lagrange multiplier increments can be trivially condensed one after another: the third and sixth row of (112) can be solved for the increments \(\Delta \varvec{\uplambda }^c\) and \(\Delta \varvec{\uplambda }^\theta \) easily, since \({{\mathbf {\mathsf{{D}}}}}_{\{u,T\}}\) are diagonal. After inserting those values for the Lagrange multiplier increments in the other lines, the remaining linear system to be solved consists of displacement and temperature degrees of freedom only:

(113)

Having solved this reduced system, the thermal and contact Lagrange multipliers can be recovered using the third and sixth row of (112), respectively, and the plastic deformation increment at every quadrature point can be computed via (63).

Numerical results

In the following, several numerical examples are presented to demonstrate the wide range of applications covered with the methods presented above. First, an anisotropic thermo-elasto-plastic necking problem without contact is analyzed. Next, two examples validate the ability of the thermo-mechanical contact algorithm to correctly model heat conductance and frictional dissipation at the interface by comparing the results to analytical and numerical reference solutions. After that, energy conservation is investigated in a dynamic thermo-mechanical contact setting. Whenever material and geometrical parameters do not model a real-world example, but only a numerical test-case, units of measurement will be omitted. Finally, a fully coupled thermo-elasto-plastic contact simulation concludes this section. All presented algorithms have been implemented in our parallel in-house multiphysics research code BACI [56].

Thermally triggered necking of an anisotropic circular bar

Necking problems are commonly analyzed with numerical methods for finite strain plasticity both in the isothermal case, see e.g. [24, 31], and in the thermo-mechanical case, see e.g. [13, 16, 17, 29]. In this contribution we want to go one step further and extend the necking problem to anisotropic thermo-plasticity. A tensile specimen with a radius of \(6.413\,\mathrm {mm}\) and a length of \(53.334\,\mathrm {mm}\) is stretched by \(26.25\%\). In the isothermal setting, this leads to a bifurcation problem, which is commonly avoided by introducing a geometric imperfection in form of an initial radius decreasing linearly towards the middle of the specimen. In the fully coupled thermomechanical setting, necking can be triggered by an inhomogeneous temperature distribution caused by convective heat transfer on the entire boundary of the specimen. The normal spatial heat flux is thereby defined as \(\mathbf {q}\cdot \mathbf {n} = -h_c (T-T_\infty )\), where \(h_c\) denotes the coefficient of convection, T the temperature at the surface, and \(T_\infty \) the temperature of the surrounding medium. Plastic anisotropy is introduced by reducing the normal yield stress in one transversal direction (direction of point \(\mathsf {A}\) in Fig. 1a) by \(17.5\%\) (see \(\alpha _1\) in Table 1). The employed elastic free energies in (33) and plastic potential in (18) are given as

All material parameters are summarized in Table 1. Exploiting the symmetries of the specimen, only one eighth of the model is discretized with 2250 first order hexahedral elements, see Fig. 1a, with appropriate symmetry conditions on the displacements and temperatures as well as a gradual mesh refinement towards the expected necking zone. Volumetric locking due to the volume-preserving plastic deformation is avoided by employing enhanced assumed strain (EAS) elements with nine additional strain modes, see e.g. [57]. When the monolithic solution algorithm presented above is combined with enhanced strain elements, coupling effects of the additional strain modes with all other fields have to be considered carefully. The modes appear as additional (element-wise discontinuous) unknowns in the system (60)–(62), having common coupling terms with the discrete displacements \({{\mathbf {\mathsf{{d}}}}}\), temperatures \({{\mathbf {\mathsf{{T}}}}}\) and plastic deformation increments \(\Delta {{\mathbf {\mathsf{{D}}}}}^p\). The local condensation procedure in (64)–(66) then becomes a two-stage process: first, at Gauss-point level, the plastic deformation increment is eliminated, and secondly, the additional strain modes are condensed at element level.

Figure 1b, c illustrate the final deformed stage and Fig. 2a the force-elongation curve for the isotropic and anisotropic thermo-mechanical necking problem. The isotropic case therein reproduces the results presented in [17]. In the first phase up to an elongation of approximately \(3.5\,\mathrm {mm}\) the deformation is dominated by homogeneous plastic deformation in longitudinal direction, such that the reaction force is dominated by plastic hardening, and the influence of anisotropy is very low. Once the necking is initiated, the plastic deformation is anisotropic in the transversal plane of the specimen, thus resulting in an anisotropic deformation pattern, see Figs. 1c and 2b. The anisotropy in the temperature distribution in Fig. 2c is less pronounced and the temperatures in points \(\mathsf {A}\) and \(\mathsf {B}\) follow the temperature evolution in the isotropic case.

Stationary heat conduction

In the next two short examples we demonstrate the accurate representation of thermal effects at the contact interface. We start with the pressure dependent heat conductivity over a non-matching contact interface, i.e. the first term in (76). This setting has already been studied in [2, 3] for node-to-segment contact formulations with a matching interface discretization and a similar setting in [4]. Contact between two elastic blocks is analyzed, where the lower surface of the lower block is supported and kept at a fixed temperature of 20, while the upper surface of the upper block is subjected to an increasing Neumann load and has a fixed temperature of 40, see Fig. 3a. Both blocks are modeled with a Saint-Venant–Kirchhoff material with Young’s modulus \(E=4000\), Poisson’s ratio \(\nu =0\), a heat conductance of 52 and no thermal expansion. For this setup, there exists an analytical solution for the steady state [3]. Figure 3b shows the resulting linear temperature distribution within each block. As can be expected for mortar methods, the contact patch test is passed to machine precision, i.e. a spatially constant contact pressure can be transmitted exactly. Moreover, Figure 3c compares the interface temperatures depending on the normal contact pressure with the analytical solution, which is recovered perfectly. Notably, the results are independent of the choice of slave and master side and do not require a matching interface discretization as used in previous studies [2, 3].

Frictionless contact with a rigid obstacle: convergence study

To further study the accuracy of the presented method, spatial convergence is investigated in this example. Since the consideration of (thermo-) plasticity does not alter the discretization compared to well-studied return mapping algorithms, but only the nonlinear solution procedure, we restrict ourselves to the case of thermo-elasticity. We study the two-dimensional problem of a rectangular thermo-elastic body \({{\mathbf {X}}}\in [-1,1]\times [0,1]\) with a rigid circular obstacle with an (outer) radius of 1.25, see Fig. 4a. Temperatures are fixed at the lower boundary of the rectangular block at \({{\mathbf {X}}}\in [-1,1]\times 0\) to \(\hat{T}=0\) and the obstacle has the fixed temperature \(\hat{T}=1\). The material is again modeled using the strain energy functions (114) and (115) with parameters \(\kappa =\nicefrac {5}{9}\) and \(G=\nicefrac {5}{12}\) (corresponding to Young’s modulus \(E=1\) and Poisson’s ratio \(\nu =0.2\)) under plane strain condition. The coefficient of thermal expansion is set to \(\alpha _T=0.01\) and heat conductivity to \(k_0=1\). Initially, both bodies are in stress-free contact; then the obstacle is moved vertically by 0.3 quasi-statically, see Fig. 4b for the deformed configuration and the resulting temperature distribution. We perform uniform mesh refinement by splitting one quadrilateral element into 4 new ones starting with \(2\times 1\) first order bi-linear elements on the coarsest level and \(256\times 128\) elements on the finest level. Since no analytical solution exists for finite deformation thermo-mechanical contact problems, the errors in displacement and temperature approximation in \(H^1(\Omega )\) (semi-) norms

are calculated with regard to a numerical reference solution \({{\mathbf {u}}}^\text {ref}\) and \(T^\text {ref}\) obtained by a mesh of \(1024\times 512\) elements. In Fig. 4c, we can observe the expected optimal convergence order of \(\mathcal {O}(h)\) for both displacement and temperature solution. Higher order approximations are not considered here, since the convergence orders are then no longer determined by the discretization itself but rather by the regularity of the solution. This holds for both contact problems [58] as well as problems of elasto-plasticity [59].

Frictional heating of a rotating ring

The second short validation example considers the effects of frictional heating. Therefore we choose an example similar to the one in e.g. [42]. A block (dimensions \(100\times 25\times 10\)) is pushed on a rotating ring (outer radius \(R_a=100\), inner radius \(R_i=75\), thickness \(z=10\)) via a constant vertical Neumann load \(F_n=150\) distributed over the top surface. Moreover, this top surface is fixed in horizontal direction and kept planar at all time. The ring is rotating at different angular velocities \(\omega \), which are applied as a Dirichlet boundary condition at the inner surface. At the contact between the block and the ring, we take the block as the slave side and set the temperature dependent friction coefficient \(\mu \) according to (73) with the parameters \(\mu _0 = 0.2\), \(T_0=293\) and \(T_d=493\). Additionally, we set the heat transfer parameters \(\bar{\gamma }^{(1)}=0\) and \(\bar{\gamma }^{(2)}=1\), which gives \(\beta _c=\delta _c=0\), meaning that there is no heat flux across the contact surface and the entire work of friction is converted to heat in the master body (i.e. the ring). By doing so, and by assuming an instantaneous heat conduction (\(k_0\rightarrow \infty \)), or equivalently low rotational speeds \(\omega \rightarrow 0\), an analytical solution for the temperature, or equivalently the thermal energy \(E_\text {thr}\), of the ring can be derived:

Frictional heating of a rotating ring. Temperature distribution after one full rotation (\(\alpha =2\pi \)) at different angular frequencies

where \(\alpha (t)\) is the rotation angle of the ring over time. To keep the focus on the thermo-mechanical contact and avoid potential thermo-elastic dissipation effects, we do not account for thermal expansion in this case, i.e. \(\alpha _T=0\) and assume the structural response to be quasi-static. Both bodies are modeled with a Saint-Venant–Kirchhoff material with Young’s moduli \(E_\text {block}=2\) and \(E_\text {ring}=10\) and Poisson’s ratios \(\nu _\text {block}=\nu _\text {ring}=0.25\). In addition, we assume \(k_0=6\) and \(C_v=10^{-3}\) for both bodies. Figure 5 displays the temperature distribution in the ring for different angular frequencies. Clearly, the higher the angular frequencies are, the more inhomogeneous the temperature distribution becomes; for the lowest shown frequency of \(\omega =10^{-2}\) an almost constant temperature across the entire ring is obtained. In that last case, the resulting change in thermal energy also agrees with the analytical solution (118) denoted as \(\mu (T)\) in Fig. 6, whereas the higher frequencies show a lower increase in thermal energy. At the onset of rotation, all curves start with the same slope, which corresponds to the analytical solution for a constant coefficient of friction (denoted as \(\mu =\mathrm {const}\)). In the following, the apparent coefficient of friction drops due to the increasing temperature, thus reducing the slope in the energy gain. Higher angular velocities result in locally higher temperatures in the contact zone and therefore a lower friction coefficient, consequently reducing the increase in thermal energy. Ultimately, all energy curves in Fig. 6 saturate at the black dashed line, which corresponds to a homogeneous temperature \(T_d\) in the entire ring. Approaching this temperature, the friction coefficient tends to zero, thus precluding any further thermal energy being introduced via work of friction at the interface.

Fig. 6

Frictional heating of a rotating ring. Change of thermal energy over rotation angle at different angular frequencies

Bouncing ball

The next example considers contact dynamics and especially matters of energy conservation in the discrete system. We consider a bouncing hollow sphere (outer radius 4, inner radius 1) between two rigid plates (dimensions \(50\times 1\times 10\), distance 14). The ball is modeled with a Saint-Venant–Kirchhoff material with Young’s modulus \(E=12.5\) and Poisson’s ratio \(\nu =0.2\). The heat capacity and conductivity of the ball and the two plates are set to \(C_{v,\text {ball}}=0.1\), \(C_{v,\text {plate}}=1\), \(k_{0,\text {ball}}=1\) and \(k_{0,\text {plate}}=0.1\), respectively. Initially, the ball is given an velocity of \(\sqrt{2}\) in a 45 degree angle towards the lower plate and a superimposed spin around an inclined axis. Initial temperatures are 1 in the ball and 0 in the plates. We set the coefficient of thermal expansion \(\alpha _T=0\), such that no transfer from mechanical to thermal energy is possible and vice-versa. The ball is discretized with 3456 second order NURBS elements and the plates with 459 elements each. As outlined in Remark 1, no special treatment of isogeometric thermo-mechanical contact is necessary compared to classical finite elements based on Lagrange polynomials, but the only difference lies in the ansatz functions employed in (43).

Fig. 7

Bouncing ball—temperature distribution and deformation state for \(\rho _{u,\infty }=\rho _{T,\infty }=0.9\) at different time steps

Figure 7 shows the deformation and temperature distribution at different time steps. During the five contact events, heat is transferred from the ball to the plates and the ball cools down to a final temperature of approximately 0.25, meaning that \(75\%\) of the thermal energy in the system is transmitted through the contact interface. Owing to the relatively coarse discretization of the plates in combination with their low thermal conductivity, slight oscillations in the temperature distribution occur in the surroundings of the contact zone, see Fig. 7b. This is due to the inability of the coarse discretization to correctly reproduce the steep gradients in the temperature field and the effect vanishes for finer discretizations. Figure 8 shows the relative change in mechanical and thermal energy over time for different spectral radii of the generalized-\(\alpha \) time integration schemes. As outlined in Remark 2, exact algorithmic energy conservation cannot be expected with the employed time integrator, which becomes obvious in the unstable behavior for \(\rho _\infty =1\), i.e. a gain of about \(8\%\) in mechanical energy after the five contact events (see Fig. 8a). Smaller values of \(\rho _\infty \) yield an energy dissipative, yet stable behavior. The thermal energy in Fig. 8b, on the other hand, is conserved to a very high accuracy, especially considering the fact that a significant amount is transmitted through the contact interface. We again want to point out that more sophisticated time integrators for contact mechanics, thermo-elasticity and elasto-plasticity are available in the literature (see references in Remark 2) but go far beyond the scope of this contribution.

Squeezed elasto-plastic tube—deformed configurations at different times including accumulated plastic strain and temperature distribution (results of an eighth model reflected for vizualization)

Squeezed elasto-plastic tube

Finally, we present a fully coupled thermo-elasto-plastic contact example to demonstrate the robustness and efficiency of the developed algorithm. Similar to the example in [20] and originally inspired by [21], a squeezed metal tube with an inner and outer radius of 4 and \(5\,\mathrm {cm}\), respectively, and a length of \(40\,\mathrm {cm}\) is analyzed. In the middle of the tube it is squeezed by two rigid cylindrical tools with an inner and outer of radius 4.5 and \(5\,\mathrm {cm}\), respectively, and a length of \(16\,\mathrm {cm}\), see Fig. 9. The material properties are again the ones given in Table 1, with plastic isotropy, i.e. \(y_{11}=y_0\). Between the tools and the tube, frictional contact with a temperature dependent friction coefficient according to (73) is assumed with the initial coefficient of friction \(\mu _0=0.25\), the reference temperature \(T_0=293\,\mathrm {K}\) and the damage temperature \(T_d=1793\,\mathrm {K}\). The tools are initially in stress free contact and perform a vertical displacement of \(u(t)=(1-\mathrm {cos}(\frac{t}{1\,\mathrm {s}}\pi ))\cdot 1.75\,\mathrm {cm}\) over time. Figure 10 illustrates the plastic strain and temperature distribution at different times. Due to the symmetry of the problem, only one eighth of the entire model is discretized with about 20.000 elements, and the results are reflected for visualization purposes. First order hexahedral elements with an F-bar technology are used to avoid volumetric locking, see [60] for the original isothermal formulation of this element. In the early deformation stages, plastic deformation and therefore heat generation is mainly located directly beneath the contact zone (see Fig. 10a), whereas later the main plastic deformation occurs at the side of the tube, where the highest peak temperatures are reached (see Fig. 10b). After contact is released, thermal conduction tends to equilibrate the temperature inhomogeneity, see Fig. 10c. To illustrate the efficient nonlinear solution procedure using Newton’s method with a consistent linearization, Fig. 11 displays the convergence behavior of different residual contributions in the time step of maximal tool velocity (\(t=0.5\,\mathrm {s}\)). All residuals clearly exhibit a quadratic rate of convergence asymptotically, until they are at some point limited by machine precision. The residual of the NCP function for plasticity (102), for instance, has been reduced by ten orders of magnitude within the first seven iterations and is then limited by numerical accuracy. In the final iteration steps, also the other residual contributions converge rapidly as expected. Finally, we compare the coupled thermo-elasto-plastic results with an isothermal structural simulation. Figure 12 displays the total contact forces for the two cases: initially, the temperature changes are low and the contact forces for both cases practically coincide. As the temperature increases, the coupled thermo-mechanical analysis softens and terminally yields contact forces about \(4\%\) lower than in the isothermal case, with a peak temperature change of \(36\,\mathrm {K}\). The presented monolithic scheme for thermo-elasto-plastic contact solves this problem within only 100 time steps with a constant step size of \(\Delta t=0.02\,\mathrm {s}\), of which the first 55 steps involve contact. Each of those time steps is solved with a standard Newton–Raphson scheme, which required an average of 9.9 iterations to converge to machine precision. The linearized system (113) is thereby solved using a preconditioned GMRES solver with a block Gauss–Seidel iteration to decouple the structural and thermal problem. Both subproblems then use incomplete factorization preconditioners [61]. However, more efficient linear solvers using algebraic multigrid (AMG) methods may be constructed based on [62] in combination with specialized AMG preconditioners for contact problems [63]. Such solvers are beyond the scope of this work, but will be a topic of future research.

Fig. 11

Squeezed elasto-plastic tube—convergence of different residuals in Newton’s method for \(t=0.5\,\mathrm {s}\). “coupled” shows the residual of (113), “structure” the residual (100), “thermo” the residual (101), “str. LM” the residual of (103) and (104), “thr. LM” the residual (105) and “plast.” the residual (102)

Conclusions

In this article, a fully coupled monolithic formulation for thermo-elasto-plasticity and thermo-elasto-plastic frictional contact problems has been derived. Following one common theme, the newly developed algorithm deals with all arising inequality constraints of contact, friction and plasticity by using nonlinear complementarity functions. The algorithm for solving the thermo-elasto-plastic problem in the bulk structure is based on our previous developments in [20], where the plastic constraints are treated via nonlinear complementarity functions at every quadrature point.

In the presented extension to thermo-plasticity, the emerging thermo-mechanical coupling terms within a monolithic solution framework are consistently linearized. By condensation of the plastic variables at quadrature level, the method has the same computational efficiency as the radial return mapping, while, at the same time, being potentially more robust, since the plastic inequality constraints only need to be satisfied at convergence of the global Newton-Raphson method, but not at every step as in the return mapping. Moreover, this general concept can be used to develop other types of complementarity functions, which may result in even more robust and efficient computational schemes.

Concerning thermo-mechanical frictional contact at finite deformations, a new discretization approach based on dual mortar finite element methods has been derived. The model includes a pressure dependent heat conduction across the contact interface as well as frictional work, which is consistently converted to heat. The use of dual basis functions allows for an easy condensation of the additional Lagrange multiplier degrees of freedom, such that the resulting system is no longer of saddle-point-type as many other contact formula are. Since the final global system to be solved only involves displacement and temperature degrees of freedom, this renders the algorithm very efficient with a quadratic rate of convergence in Newton’s method, yet at the same time very accurate as the obtained spatial convergence results underline. Finally, robust compatibility of the presented methods has been demonstrated in a fully coupled thermo-elasto-plastic frictional contact setup.

Declarations

Author's contributions

AS developed the idea, conducted numerical experiments and wrote draft. WAW and AP fine-tuned the research idea, suggested numerical experiments and revised the paper. All authors read and approved the final manuscript.

Acknowledgements

None.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

Not applicable.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Funding

Funding was provided by the German Research Foundation (Deutsche Forschungsgemeinschaft—DFG) in the framework of the Priority Programme “Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretisation Methods, Mechanical and Mathematical Analysis” (SPP 1748) under projects PO1883/1-1 and WA1521/15-1.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.