Numerical analysis of a nonmonotone dynamic contact problem of a non-clamped piezoelectric viscoelastic body

We consider a contact process between a body and a foundation. The body is assumed to be viscoelastic and piezoelectric and the contact is dynamic. Unlike many related papers, the body is assumed to be non-clamped. The contact conditions has a form of inclusions involving the Clarke subdifferential of locally Lipschitz functionals and they have nonmonotone character. The problem in its weak formulation has a form of two coupled Clarke subdifferential inclusions, from which the first one is dynamic and the second one is stationary. The main goal of the paper is numerical analysis of the studied problem. The corresponding numerical scheme is based on the spatial and temporal discretization. Furthermore, the spatial discretization is based on the first order finite element method, while the temporal discretization is based on the backward Euler scheme. We show that under suitable regularity conditions the error between the exact solution and the approximate one is estimated in an optimal way, namely it depends linearly upon the parameters of discretization.


1.
Introduction. Contact processes appear in many fields of industry and everyday life. Hence it is desirable to study mathematical models concerning various properties and parameters of physical bodies remaining in contact. In particular, the piezoelectric effects play a special role in such models. It is enough to mention that the use of the smartphone screen is based on piezoelectric contact and most of contemporary human activities would be impossible without this phenomenon. There are numerous publications concerning mathematical analysis of contact processes with piezoelectric effect, that essentially deal with existence and uniqueness of solution for the studied problems. We refer, for example, to [19][20][21][22]26] for some results in this field.
Apart from the theoretical aspect of mathematical modelling, the numerical analysis plays also a significant role. Typically, it includes approximate schemes based on temporal and spatial discretization, estimate of the error between the approximate and the exact solution, and computer simulations. This kind of approach carried out for various contact problems of mechanics with or without piezoelectric effect is presented in [1-7, 9, 10, 14, 17, 18], for instance.
In this paper we consider a viscoelastic body with piezoelectric properties. The body is in a contact with a foundation. We assume that the contact process is dynamic and the contact conditions are modelled by possibly nonmonotone multivalued relations involving the Clarke subgradient of locally Lipschitz functionals. A similar model was studied in [20], where the existence and uniqueness of a weak solution have been proved. However, the body considered in [20] is assumed to be clamped which prevents it from moving away. In other words, the Dirichlet boundary condition is assumed at the part of body's surface. This assumption, which is usually used in many papers, allows to simplify mathematical analysis of the problem by using the Korn inequality. Thanks to it, the space of test functions used in the variational formulation, may be equipped with a norm which involves only the strain tensor of the displacement function. In contrast, we deal with a non-clamped body and we can not use the Dirichlet boundary condition. Hence, in our case, the corresponding norm involves not only the deformation but also the displacement function itself. In a consequence, the numerical analysis of our problem requires more sophisticated technique than the clamped case, as we have more terms to estimate.
In this paper we do not prove existence result for the studied problem, as we concentrate on its numerical analysis. For the approximation of the original problem we use a fully discrete numerical scheme in which the temporal discretization is based on the backward (implicit) Euler scheme, while the spatial discretization is based on the first order finite element method. We provide an optimal bound for an error between the solution of original problem and the approximate one with respect to the discretization parameters.
The presented result generalizes the one obtained in [3], where the nonmonotone contact process of a clamped viscoelastic body has been studied and the numerical analysis has been carried out. The generalization of [3] is twofold. First, as mentioned before, the body under consideration is non-clamped, and second, an additional piezoelectric property of the body is considered. For other results concerning numerical analysis of a non-clamped body, we refer, for example, to [9]. The rest of the paper is structured as follows. In Section 2 we provide notation and preliminary material to be used later. In Section 3 we present a mathematical model and describe its physical background. Next, we list assumptions on its data and provide its variational formulation. In Section 4 we construct a fully discrete scheme and carry out the numerical analysis. The main result of this section, Theorem 4.2, provides an optimal error estimate for the numerical solution.
2. Notation and preliminaries. In this section we present the notation and some preliminary material. For further details, we refer the reader to [13,15,16,25].
We denote by S d the space of second order symmetric tensors on R d (d ≤ 3 in applications). The corresponding inner products and the norms on R d and S d , respectively, are given by Here and below, the indices i and j run between 1 and d, and the summation convention over repeated indices is adopted.
Let Ω ⊂ R d be a bounded domain with a Lipschitz boundary Γ and ν denote the unit outward normal on Γ. We introduce the following function spaces.
Here ε : V → Q, Div : Q 1 → H and div : H → L 2 (Ω) are the deformation and divergence operators, defined by respectively, where the index following a comma indicates the partial derivative with respect to the corresponding component of the independent variable. The spaces H, Q, V , Q 1 and W are real Hilbert spaces endowed with the canonical inner products given by The associated norms on these spaces are denoted by · H , · Q , · V , · Q1 and · W , respectively. Moreover, we use notation We recall that | · | V is a seminorm on the space V . Letγ : V → L 2 (∂Ω) be the trace map. For every element v ∈ V we use the same symbol v to denote the traceγv of v on Γ, and we denote by v ν and v τ the normal and tangential components of v on the boundary Γ given by We denote by σ ν and σ τ the normal and tangential traces of σ, Now we recall two versions of the Green formula (cf. [23], (2.6) and (2.7), respectively) needed to obtain a variational formulation of the problem under consideration.
Next we recall the definitions of the classical (one-sided) directional derivative and its generalization in the sense of Clarke. Let X be a Banach space and X * be its dual. For a function ϕ : X → R the directional derivative of ϕ at x ∈ X in the direction v ∈ X is defined by whenever this limit exists. The Clarke generalized directional derivative of a locally Lipschitz function ϕ : X → R at the point x ∈ X in the direction v ∈ X, denoted by ϕ 0 (x; v), is defined by The Clarke subgradient of ϕ at x, denoted by ∂ϕ(x), is a subset of X * given by A comprehensive introduction to the theory of the Clarke subgradient can be found in classical handbooks [12,24]. We also refer to [23] for more recent results in this field. We will need the following discrete Gronwall inequality ( [15,Chap. 7]). We complete this section with the following result which will be used in next sections.
The proof of Lemma 2.2 follows the lines of the proof of Lemma 30 in [8].
3. Mechanical problem and its variational formulation. We start with a description of the physical problem. A viscoelastic body made of a viscoelastic, piezoelectric material occupies an open, bounded, connected set Ω ⊂ R d (d = 2, 3) with a Lipschitz boundary Γ. We consider a partition of the boundary into two disjoint parts Γ N and Γ C with Γ N and Γ C being relatively open. On the other hand, the part Γ N is divided into two parts Γ a and Γ b , such that meas (Γ a ) > 0. A volume force of density f 0 acts in Ω and a surface traction of density f N is applied on Γ N . The body has volume electric charge of density q 0 . We also assume that the electrical potential vanishes on Γ a and a surface charge of density q b is prescribed on Γ b . The body is in contact with an obstacle on Γ C . The obstacle is assumed to be electrically conductive and its potential is maintained at ϕ F . The contact is frictional and there may be electrical charges on the contact surface. Let [0, T ] be a time interval of interest, T > 0. For any time dependent function v, we denote byv andv its first and second order time derivatives, respectively.
The classical formulation of the mechanical problem is the following.
Now we provide a short description of conditions (4)- (14). For simplicity, we skip the dependence of given functions and operators upon spatial and time variables. Equations (4) and (5) represent the electro-viscoelastic constitutive law, in which ε(u) denotes the linearized strain tensor, A is the viscosity operator, assumed to be nonlinear, B is the elasticity tensor, P is the third-order piezoelectric tensor, P is its transpose, β denotes the electric permittivity tensor and E(ϕ) is the electric field. We recall that E(ϕ) = −∇ϕ. The tensors P = (p ijk ) and P satisfy the equality and the components of tensor P are given by (P ) ijk = p kij . Condition (6) is the equation of motion in which the density of mass has been taken to equal one and (7) represents the balance equation for the electric-displacement field. Equation (8) is the traction boundary conditions, which expresses the fact that the density of the boundary force is prescribed on Γ N . Moreover, (9) and (10) represent the electric boundary conditions. They model the fact that the electric potential vanishes on Γ a , while the electric charges are prescribed on Γ b . Now we pass to the description of contact conditions. In particular, (11) and (12) represent the dependence of normal and tangential stresses upon normal and tangential velocities, respectively. Moreover, (13) represents electrical conductivity condition. The symbol ∂ used in relations (11)-(13) stands for the Clarke subgradient of a locally Lipschitz potential.
In general the potentials j ν , j τ and j may be nonconvex and hence their Clarke subgradients may be nonmonotone. In a consequence, we deal with nonmonotone contact conditions. Finally, (14) represents the initial conditions. In the study of Problem P M we need the following assumptions on its data.

H(A): The viscosity operator
.
We remark that the hypotheses H(j ν )(b), H(j τ )(b) and H(j)(b) are used in order to guarantee that the subgradients ∂j ν , ∂j τ , ∂j are well defined. But in fact, due to the Lebourg mean value theorem, the hypotheses H(j ν )(c), H(j τ )(c) and H(j)(c) imply that the functionals j ν , j τ and j are globally Lipschitz, which is stronger Finally we impose the following assumptions on the constants m β and m j introduced in assumptions H(β) and H(j), and the trace operatorγ introduced in Section 2.
Let us define the following constants.
In this paper we do not deal with existence of solution to Problem P V . However, we remark that Theorem 3 of [20] provides existence of unique solution for a problem analogous to Problem P V with the regularity u ∈ L 2 (0, T ; V ),u ∈ L 2 (0, T ; V ), u ∈ L 2 (0, T ; V * ) and ϕ ∈ L 2 (0, T ; Φ). The mechanical process studied in [20] is assumed to be clamped, namely there is a part of boundary Γ D ⊂ ∂Ω of positive measure such that u = 0 on Γ D . Moreover, an additional smallness assumption for the constants m ν and m τ is used there. Both these restrictions make the problem easier to analyse. According to the authors knowledge, the hypotheses listed in this section are sufficient to obtain the analogous existence and uniqueness result for the non-clamped problem presented here, even without any additional smallness assumption. This is going to be a subject of a future dissertation. 4. Discretization and error estimates. In this section we consider an approximate problem corresponding to Problem P V based on the finite element discretization and provide an error estimate between the solutions of both problems. For simplicity we assume that Assume that Ω is a polygonal domain in R 2 or a polyhedral domain in R 3 . Then, each part of the boundary Γ N and Γ C is a sum of segments (in the case Ω ⊂ R 2 ) or polygons (if Ω ⊂ R 3 ), having mutually disjoint interiors, i.e., Let {T h } be a family of regular triangulations of the setΩ of diameter h (cf. [11], Chapter 3, § 3.1) which introduces a division ofΩ into triangles/tetrahedrons compatible with the division of the boundary ∂Ω into parts Γ j,i , 1 ≤ i ≤ N j , j ∈ {N, C}, in that sense that if any edge of a triangle/any face of tetrahedron that belongs to T h has an intersection of positive boundary measure with any of sets Γ j,i , then the whole edge/face is contained in Γ j,i . For an element T ∈ T h let P 1 (T ) denote the space of polynomials of order ≤ 1 on T . We introduce the following spaces of continuous, piecewise affine functions on elements T ∈ T h .
Let Π h denote the finite element interpolant (see [11], (2.3.29)). For the sake of simplicity, we denote by the same symbol Π h both operators Π h : V → V h and Π h : Φ → Φ h . In what follows, we recall the classical approximation properties of Π h in the spaces V and Φ. Let u ∈ V ∩ H 2 (Ω; R d ) and ϕ ∈ Φ ∩ H 2 (Ω). Then, we have with a constant c independent of h. Additionally, we define a uniform partition of [0, T ] denoted by 0 = t 0 < t 1 < . . . < t N = T , where N ∈ N, t j = jk for j = 1, ..., N and k = T N is the time step. For a continuous function h we use the notation g n = g(t n ). For a sequence {z n } N n=0 we denote by δz n = (z n − z n−1 )/k for n = 1, . . . , N its backward divided difference.

KRZYSZTOF BARTOSZ
Let u hk 0 = Π h u 0 and u hk 1 = Π h u 1 . The fully discrete approximation of Problem P V is the following.
Let the triple u, w and ϕ be a solution of Problem P V and the triple {u hk n } N n=1 , {w hk n } N n=1 and {ϕ hk n } N n=1 be a solution of Problem P hk V . Let us introduce the notation In what follows, we will denote by c a generic positive constant independent on the discretization parameters h and k which can differ from place to place. In contrast, we will denote by c( ) a constant that depends on a parameter . We are now in a position to formulate the theorem concerning the estimate of an error between the solutions of both problems.
there exists a positive constant c independent of j such that Proof. We deal with the proof of (40) first. We test equation (21) at a fixed point t j , j = 0, 1, ..., N by an element ψ h ∈ Φ h . We also test equation (36) with n = j by the same element and subtract it from (21). Then, we obtain Since (42) holds for an arbitrary element ψ h ∈ Φ h , in particular, we have In what follows, we estimate the left hand side of (43) from below and its right hand side from above using assumptions H(β), H(P) and H(j). In particular, we get for Using condition (16), we can choose such > 0 that m β −m j γ 2 −2 is positive. Moreover, we replace ψ h by ψ h j and obtain (40). Now we pass to the proof of (41). To this end we consider equation (17) at the point t j , j = 1, ..., N and test it by an element v h ∈ V h . We test equation (32) with n = j by the same element and subtract it from (17). Hence, we get Note that (50) holds for arbitrary v h ∈ V h . Hence, in particular, we have After some reformulation, we obtain In what follows, we estimate the left hand side of (51) from below and its right hand side from above using assumptions H(A), H(B) H(P), H(j ν ) and H(j τ ). In particular for the left hand side, we get for any > 0 We remark that the derivation of (56) and (57) is based on Lemma 2.2. As for the right hand side of (51), we get Applying (52)-(63) to (51) we get 1 2k We add the term (m A − (3 + m ν + m τ )) w j − w hk j 2 H to both sides of (64). Next, we substitute (40) to (64) and replace v h by v h j . Then, we obtain 1 2k . Now we sum up inequalities (65) for j ∈ {1, ..., n}, n = 1, ..., N and multiply the resulting inequality by 2k. Then, we obtain Let < m A (3 + m ν + m τ ) −1 . Then, it follows from (66) that Furthermore, using (34), we can write where I j is defined by (38). Hence, we have Next, we calculate Substituting (69)-(70) to (67), we get We introduce the notation e n = w n − w hk n 2 H + k n j=1 w j − w hk j and g n = w 0 − w hk Then (71) is equivalent with Hence, applying Lemma 2.1 to (71), we obtain (41), which completes the proof.
Theorem 4.1 is a starting point to obtain the optimal error estimates for the numerical scheme. However, we need the following a-priori assumptions on the regularity of solution to Problem P V .
where the constant c > 0 does not depend on h and k.
Proof. Recall that u hk 0 = Π h u 0 and w hk 0 = Π h w 0 . Moreover, we take v h j = Π h w j and ψ h j = Π h ϕ j in (41). By approximation properties (26)-(31), we estimate each term of the right hand side of (41).
For the remaining terms, we have Hence, it follows that