Modelling and simulation of a dynamic contact problem in thermo-piezoelectricity

In this work, we numerically study a dynamic frictional contact problem between a thermo-piezoelectric body and a conductive foundation. The linear thermo-electro-elastic constitutive law is employed to model the thermo-piezoelectric material. The contact is modelled by the Signorini condition and the friction by the Coulomb law. A frictional heat generation and heat transfer across the contact surface are assumed. The heat exchange coefficient is assumed to depend on contact pressure. Hybrid formulation is introduced, it is a coupled system for the displacement field, the electric potential, the temperature and two Lagrange multipliers. The discrete scheme of the coupled system is introduced based on a finite element method to approximate the spatial variable and an Euler scheme to discretize the time derivate. The thermo-mechanical contact is treated by using an augmented Lagrangian approach. A solution algorithm is discussed and implemented. Numerical simulation results are reported, illustrating the mechanical behavior related to the contact condition.


Introduction
T he effective conversion of the electrical energy into mechanical energy and vice versa has led the piezoelectric materials to important applications in many engineering structures such as sensors, actuators, intelligent structures, etc. Thermal effects, such as temperature induced deformation and the pyroelectric effect, are especially important for many smart ceramic materials. Thus, a coupling of thermo-electro-mechanical fields is needed to be taken into account if a temperature load is considered in a piezoelectric solid. Models taking into account thermal and piezoelectric effects can be found in [1][2][3].
Thermal effects in contact processes affect the composition and stiffness of the contacting surfaces, and cause thermal stresses in the contacting bodies. When extending contact problems to thermomechanics, additional thermal effects need to be accounted for: heat conduction appears through the contact interface and frictional work is converted to heat. Recent models of frictional contact problems involving thermo-piezoelectric materials can be found in [4][5][6] and the references therein. There, besides the rigorous construction of various mathematical models of contact for thermo-piezoelectric materials, the unique weak solvability of these models was proved, by using arguments of variational and hemivariational inequalities. Numerical analysis of the problems, including numerical simulations, can be found in [7,8]. In [7] the process was assumed to be static, the contact was described with Signorini condition and Tresca's law of dry friction, and a regularized electrical and thermal conductivity conditions and, in [8], the process was assumed to be quasistatic, the contact was assumed to be bilateral, in which contact is always maintained; and associated to the Tresca friction law, and to the heat exchange condition. There, discrete schemes to approximate the problems were considered and implemented in a numerical code, and numerical simulations were provided.
The present paper represents a continuation of [8] and it deals with a mathematical model which describes the frictional contact between a thermo-piezoelectric body and a thermally conductive foundation. We use both the thermo-electro-elastic constitutive law used in [8] but unlike [8], we assume here that the process is dynamic and the foundation is completely rigid and we model the contact with the Signorini condition with Coulomb's law of dry friction. This condition is other physical setting (see, e.g., [9,10]). Also, note that, unlike [8]; here the model includes frictional heat generation and the heat exchange condition, in which the heat exchange coefficient is not a constant but a function of the contact pressure (see [11,12]). The other trait of novelty of the present paper consists in the fact that here we deal with the numerical approach of the problem and provide numerical simulations. The corresponding numerical scheme is based on the spatial and temporal discretization. Furthermore, the spatial discretization is based on the finite element method, while the temporal discretization is based on the Euler scheme. Then, the scheme was utilized as a basis of a numerical code for the problem, in which we develop a specific contact element. We need this element in order to take into account the coupling of the mechanical and thermal unknows on the contact boundary condition. By using the code, simulation results on numerical example are presented.
The rest of paper is structured as follows. In Section 2, we describe our model. Section 3 introduces the variational formulation of the problem, and a fully discrete variational approximation by considering a hybrid formulation. The numerical algorithm used for solving the discrete problem is described in Section 4, where some numerical simulations are also presented to highlight the performance of the method and the effects of the conductivity of the foundation, as well. Finally, in Section 5, we present some conclusions and perspectives.

The model
The physical setting is the following: A piezoelectric body occupies a regular domain Ω ⊂ R d (d = 2, 3) with a smooth boundary ∂Ω = Γ. The body is submitted to the action of body forces of density f 0 , a volume electric charges of density φ 0 and a heat source of constant strength ϑ 0 . It is also submitted to mechanical, electric and thermal constraints on the boundary. To describe them, we consider a partition of Γ into three measurable parts Γ D , Γ N , Γ C , on one hand, and a partition of Γ D ∪ Γ N into two open parts Γ a and Γ b , on the other hand. We assume that meas Γ D > 0 and meas Γ a > 0. The body is clamped on Γ D , therefore, the displacement field vanishes there. Moreover, we assume that a density of traction forces, denoted by f N , acts on the boundary part Γ N . We also assume that the electrical potential vanishes on Γ a and a surface electric charge of density φ b is prescribed on Γ b . We suppose that the temperature vanishes in Γ D ∪ Γ N . In the reference configuration, the body is in contact over Γ C with a thermally conductive foundation. The contact is modelled with a Signorini's conditions and a version of Coulomb's law dry friction. Here, we study the evolution of the state of the system on a finite time interval [0, T], with T > 0. Here and everywhere in this paper i, j, k, l run from 1 to d, summation over repeated indices is implied and the index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable, i.e., We denote by u = (u i ) ∈ R d the displacement vector, by σ = (σ ij ) ∈ S d the stress tensor, by "(u) = (ε ij (u)) ∈ S d the linearized strain tensor, i.e., ε ij (u) = (u i,j + u j,i )/2, by E(ϕ) = −∇ϕ = −(ϕ ,i ) ∈ R d the electric vector field, where ϕ ∈ R is the electric potential, and by θ ∈ R the temperature. The notation S d stands for the space of second order symmetric tensors on R d . We also use the dot to denote the time derivative, sou = (u i ) represents the velocity vector. The body is assumed to be thermo-electro-elastic and, therefore, we use the constitutive law and the heat flux vector q = (q i ) ∈ R d is given by the Fourier law of heat conduction Here D = (D i ) is the electric displacement field and F = ( f ijkl ), E = (e ijk ), β = (β ij ), M = (m ij ), P = (p i ) and K = (k ij ) are respectively, the elasticity, piezoelectric, electric permittivity, thermal expansion, pyroelectric and thermal conductivity tensors. E T denotes the transpose of E ; also, the tensors E and E T satisfy the equality and the components of the tensor E T are given by e T ijk = e kij . Since the process is assumed dynamic, then the equation of stress equilibrium, the equation of the quasistationary electric field, and the heat conduction equation are , Here α is given as α = ρc ν /θ re f , where ρ is the mass density, c ν is the specific heat and θ re f is the reference uniform temperature of the body. Moreover, Div and ÷ represent the divergence operators for tensor and vector functions, i.e., We turn to describe the boundary conditions, so we denote by ν = (ν i ) the unit outward normal on Γ. Then, on the Γ D ∪ Γ N portion of the boundary, we impose the following conditions The boundary conditions for the electric potential can be defined in the following forms: Next, on Γ D ∪ Γ N we prescribe a Dirichlet condition for the temperature, say, We now describe the thermo-mechanical boundary conditions on the potential contact surface Γ C . We assume that the normal displacement u ν = u · ν and the normal stress σ ν = σν · ν satisfy the Signorini's contact conditions The corresponding Coulomb law of dry friction may be stated as follows: Hereu τ =u −u ν ν is the tangential velocity, σ τ = σν − σ ν ν represents the tangential force on the contact boundary and µ ≥ 0 is the coefficient of friction. Next, we describe the boundary condition for the temperature on Γ C . We assume that there is heat exchange between the surface and the foundation, which is at temperature θ f . Moreover, since the flux of heat generated by the friction traction on the contact surface is proportional to the tangential shear σ τ and to the tangential velocityu τ of the surface, we assume a boundary condition of the following form where k c (·) is the normal pressure dependent heat exchange coefficient and has to satisfy k c (0) = 0. This condition guarantees that there is no heat flux between the body and the foundation if they are not in contact. For k c (·), we employ a linear model k c (σ ν ) =k c |σ ν |, wherek c ≥ 0 is model constant, see [11]. We collect the above equations and conditions to obtain the following mathematical problem. Problem P. Find a displacement field u : Divσ Here, u 0 , v 0 , ϕ 0 and θ 0 are the prescribed initial displacement, velocity, electric potential and temperature, respectively.

A hybrid variational formulation
We now turn to the variational formulation of Problem P which is the starting point for the numerical modelling based on the finite element discretization. We denote in the sequel by " · " and · the inner product and the Euclidean norm on the spaces R d and S d . We introduce the spaces and we use the notation H = [L 2 (Ω)] d , and we introduce the spaces The spaces H, V, W, Q and H are real Hilbert spaces endowed with the canonical inner products given by We consider the trace spaces X ν = {w ν| Γ C : w ∈ V} and X τ = {w τ | Γ C : w ∈ V} equipped with their usual norms. Denote respectively by X * ν and X * τ the dual of the spaces X ν and X τ . We consider the three mappings f : Then, the hybrid variational formulation of the contact problem P obtained by multiplying the equations with the relevant test functions and performing integration by part, is as follows. Problem P V . Find a displacement field u : [0, T] −→ V, a normal stress λ ν : [0, T] −→ X * ν , a tangential stress λ τ : [0, T] −→ X * τ , an electric potential field ϕ : [0, T] −→ W and a temperature field θ : [0, T] −→ Q such that for a.e., t ∈ (0, T) (ρü(t), w) H + (F "(u(t)), "(w)) H + (E T ∇ϕ(t), "(w)) H − (Mθ(t), "(w)) H The inclusion in Equation (16) represents the Signorini contact condition (9). Here, ∂ denotes the subdifferential operator in the sense of convex analysis and I R − denotes the indicator function of the negative half-line. Recall also that, the inclusion in Equation in (17) represents the subdifferential form of Coulomb's law of dry friction (10), see [9] for details.

Numerical approximation
This section is devoted to the numerical discretization of the of problem P V . First, we consider tree finite dimensional spaces V h ⊂ V, W h ⊂ W and Q h ⊂ Q approximating the spaces V, W and Q, respectively, in which h > 0 denotes the spatial discretization parameter, and let U h = U ∩ V h . In the numerical simulations presented in the next section, V h , W h and Q h consist of continuous and piecewise affine functions, that is, where Ω is assumed to be a polygonal domain, T h denotes a finite element triangulation of Ω, and P 1 (Tr) represents the space of polynomials of global degree less or equal to one in Tr.
To discretize the time derivatives, we use a uniform partition of [0, T], denoted by 0 = t 0 < t 1 < . . . < t N = T. Let k be the time step size, k = T/N, and for a continuous function f (t) let f n = f (t n ). Finally, for a sequence {w n } N n=0 we denote by δw n = (w n − w n−1 )/k the divided differences. We now consider the spaces X h its usual norm. We also consider the discrete space of piecewise constants X * ν h ⊂ L 2 (Γ C ) and X * τ h ⊂ L 2 (Γ C ) d related to the discretization of the normal and the tangential stress, respectively.
The fully discrete approximation of Problem P V , based on the Euler scheme, is the following: h , a discrete electric potential ϕ hk = {ϕ hk n } N n=0 ⊂ W h and a discrete temperature θ hk = {θ hk n } N n=0 ⊂ Q h such that, for all n = 1, . . . , N, Here, u hk 0 , v hk 0 , ϕ hk 0 and θ hk 0 are appropriate approximation of the initial condition u 0 , v 0 , ϕ 0 , θ 0 and λ hk 0 = 0.

A solution algorithm
We now describe the numerical solution of the hybrid variational Problem P hk V . The numerical treatment of the conditions (18) and (19) is based on the augmented Lagrangian approach (see [10,13] for more details). To this end we introduce the notation λ = λ ν ν + λ τ , where λ ν = λ.ν and λ τ = λ − λ ν ν. Let N h tot be the total number of nodes and denote by α i , β i , γ i the basis functions of the spaces V h , W h and Q h , respectively, for i = 1, . . . , N h tot . Then, the expression of functions w h ∈ V h , ξ h ∈ W h and η h ∈ Q h is given by where w i , ξ i and η i represent the values of the corresponding functions w h , ξ h and η h at the i th node of T h . It can be shown that the numerical approach of Problem P hk V is governed at each time step n by a system of non-linear equations of the form R(δv n , δϕ n , δθ n , v n , u n , ϕ n , θ n , λ n ) =M(δv n ) +Ã(v n , δϕ n , δθ n ) +G(u n , ϕ n , θ n ) +F(u n , θ n , λ n ) = 0, where the functionsM,Ã,G andF are defined below. Here, the vectors δv , u i n , ϕ i n , and θ i n represent the value of the function δv hk n , δϕ hk n , δθ hk n v hk n , u hk n , ϕ hk n and θ hk n at the i th nodes of T h . λ i n denotes the value of λ hk n at the i th node of the discretized contact interface, where N h Γ C denotes the total number of nodes of i th lying on Γ C . Next, the generalized acceleration termM(a) Γ C and the generalized thermo-electro-elastic term tot denotes the acceleration term, the damping term and the thermo-electro-elastic term, given by Above, w, ξ and η represent the generalized vectors of components w i , ξ i and η i , for i = 1, · · · , N h tot , respectively. Finally, the generalized contact operatorF(u, where ∇ x represents the gradient operator with respect the variable x. Also L r denote the augmented Lagrangian functional for the contact and friction terms, where r > 0 is an augmentation parameter, P C[ρ] is the orthogonal projection on the Coulomb convex disk of constant radius ρ, and (·) − is the negative part of x ∈ R, i.e., (x) − = max(−x, 0).
Γ C the thermo-mechanical frictional contact operator defined through the relation The solution of the non-linear system (20) is based on a linear iterative method similar to that used in the Newton method, which permits to treat simultaneously the four unknowns u n , ϕ n , θ n and λ n and, for this reason, we use in what follows the notation x n = (u n , ϕ n , θ n , λ n ). This Newton algorithm can be summarized by the following iteration process ); i and n represent respectively the Newton iteration index and the time index is the contact tangent matrix; also, D u M, D u,ϕ,θ A, D u,ϕ,θ G and D u,θ,λ F denote the differentials of the functions M, A, G and F with respect to the variables u, ϕ, θ and λ. This leads us to solve the resulting linear system Note that formulation (20) has been implemented in the open-source finite element library GetFEM++ (see [14]).

Numerical results
For the numerical simulations we consider the physical setting depicted in Figure 1 Here, we use as material the thermo-piezoelectric body whose constants are taken as [2]. The following data have been used in the numerical simulations: Our interest in this example is to study the influence of the thermal conductivity of the foundation on the contact process and, to this end, we consider the problem both in the case when the foundation is insulated there are no heat flux on Γ C (i.e. q · ν = 0 on Γ C ) and in the case when it is thermally conductive. Figure    2 presents the deformed configurations for the two previously mentioned cases, at final time. We can easily note that considering a thermally conductive foundation reduce the deformations. In order to highlight the influence of the foundation temperature on the electric potential, we plot the electric potential for the two previously mentioned cases (see Figure 3). The first case illustrates the direct piezoelectric effect: the electric potential is generated because of the deformation. However, in the second case, we can easily note that considering a thermally conductive foundation increases the electric potential. In Figure 4, the Von Mises stress norm is plotted on the deformed configuration. Clearly, effects due to the influence of foundation temperature, can be observed. Both temperature is plotted in Figure 5 at final time for the value θ f = 393 K.

Conclusion
In this paper thermo-piezoelectric contact including frictional heat generation and interfacial heat transfer is numerically studied. The novelties arise in the fact that the process is dynamic, the material behavior is described by a thermo-electro-elastic constitutive law and the foundation is thermally conductive. A fully discrete scheme was used to approach the problem and a numerical algorithm which combine the augmented Lagrangian approach with the Newton method was implemented. Moreover, numerical simulations for a representative two-dimensional example were provided. These simulations describe the thermal effect, i.e. the appearance of strain and voltage in the body, due to the action of the temperature field. Also, they underline the effects of the thermal contact, i.e. heat transfer and frictional heating, on the process. Performing these  simulations, we found that the numerical solution worked well and the convergence was rapid. This work opens the way to study further models of frictional contact with a coefficient of friction depending on the slip rate or the temperature.
Conflicts of Interest: "The author declares no conflict of interest."