A Contact Model for Piezoelectric Beams

. We consider a mathematical model which describes the equilibrium of an electro-elastic beam in contact with an electrically conductive foundation. The model is constructed by coupling the beam equation with the one dimensional piezoelectricity system obtained in [13]. We state the unique weak solvability of the model as well as the continuous dependence of the weak solution with respect to the data. We also introduce a discrete scheme for which we perform the numerical analysis, including convergence and error estimates results. Finally, we present numerical simulations in the study of a test problem.


Introduction
Piezoelectric materials belong to the family of "smart materials" and are characterized by a coupling of mechanical and electrical properties.Thus, electric charges can be observed on a piezoelectric body subjected to the action of external forces and, conversely, an electric potential applied on a piezoelectric body gives rise to stresses and strains.These properties of piezoelectric materials make them suitable to be used as sensors and actuators in various industrial settings and real-world applications.For this reason, the interest in the analysis of mathematical models with piezoelectric materials is currently increasing.
The construction of appropriate models to describe the behaviour of thin deformable bodies like plates, shells and beams, represents an important topic in Solid Mechanics.By using asymptotic analysis, several classical reduced models have been mathematically justified over the years.Pioneering work for modelling The work of Á. Rodríguez-Arós was supported by the project "Modelización y simulación numérica de sólidos y fluidos en dominios con pequeñas dimensiones.Aplicaciones en estructuras, biomecánica y aguas someras", MTM2012-36452-C02-01 financed by the Spanish Ministry of Economía y Competitividad with the participation of FEDER of thin linearly isotropic piezoelectric beams was performed in [1,5,8,[13][14][15]17].Models for elastic beams in contact with a foundation have been justified in [7,9,16], based on the ideas in [12].
The present paper represents a continuation of our previous works.Here we analyse, both mathematically and numerically, a model for an elastic piezoelectric beam in contact with a deformable conductive foundation.The manuscript is structured as follows.In Section 2 we describe the physical setting together with the corresponding mathematical model.Then, we list the assumptions on the data and state our main results in the analysis of the model.In Section 3 we formulate the discrete problem by using Finite Elements and provide existence, uniqueness, convergence and error estimates results.We also provide a brief description of the corresponding numerical algorithm.Finally, in Section 4 we present numerical simulations which highlight the performances of the algorithm and describe the effects of the different parameters on the solution.

Problem Statement
We consider an elastic piezoelectric beam of length L > 0, cross-section area A, Young Modulus E and Inertia Moment I.The beam is clamped at both ends and subjected to axial and vertical external forces f and f ⊥ , respectively.As a result, it may enter in contact with a conductive deformable foundation.Based on our previous works [9,13,16], we associate to this physical setting the following mathematical model.Problem 1. Find a bending field ξ : [0, L] → IR, a stretching field u : [0, L] → IR and an electric potential q : [0, L] → IR such that Here, P and ε > 0 denote the piezoelectric coefficient and the electric permittivity coefficient, respectively, and C ⊥ = EI > 0, C = EA > 0.Moreover, p(•) : IR → IR + denotes the normal compliance function, which vanishes for negative arguments, and µ 1 > 0, µ 2 > 0 represent coefficients of the system.The dependence of the various functions on x ∈ [0, L] is not indicated explicitly and the symbol stands for the derivative with respect to this spatial variable.We now briefly describe the equations and conditions (2.1)-(2.6).First, equation (2.1) is the beam equation in normal compliance contact with a foundation, with an initial gap s.It also contains an additional term on the right hand side, which describes the electric charges from the obstacle to the beam, when the contact arises.Here we use the notation R(z) = (R(z)) + , for all z ∈ IR, where r + = max {r, 0} denotes the nonnegative part of r and R is the truncation operator given by M > 0 being a positive constant that depends on the characteristic length of the system.Equations (2.2)-(2.3)are the piezoelectric equations for elastic beams presented in [13], and (2.4)-(2.6)represent the boundary conditions, in which q 0 and q L denote the electric potentials applied on both ends of the beam.
For the analysis of Problem 1 we use the standard notation for Sobolev and Lebesgue spaces.In particular, we denote by • 0 , • 1 , • 2 and • ∞ the norms on the spaces L 2 (0, L), H 1 (0, L), H 2 (0, L) and L ∞ (0, L), respectively.Let ϕ = q − q, where q is a lifting function of q 0 and q L in H 1 (0, L) (see [13]).Note that q ∈ C 0 ([0, L]) and max Moreover, without loss of generality, we assume that q 0 , q L > 0. As a consequence, it follows that q > 0 as well.We also assume that and we denote η = ξ − s.Then, using a standard procedure, it is easy to obtain the following variational formulation of problem (2.1)-(2.6).
In the study of Problem 2 we assume the following hypotheses. ) Moreover, we assume that the normal compliance function p satisfies for all r ∈ R, (e) The mapping p(•, r) : x → p(x, r) vanishes for all r ≤ 0, (2.14) and, in addition, Next, we consider the space X(0, L) = H 2 0 (0, L) × H 1 0 (0, L) × H 1 0 (0, L), which is a real Hilbert space with the canonical inner product, denoted by (•, •) X(0,L) .We then define the operator A : X(0, L) → X(0, L) and the element F ∈ X(0, L) by equalities (2.17) Then, an equivalent formulation of Problem 2 is as follows.
Our main result in the study of Problem 3 is the following.
Theorem 1.Under the assumptions (2.7), (2.11)-(2.15),there exists a unique solution x = (η, ϕ, u) ∈ X(0, L) to Problem 3.Moreover, if represents the solution of Problem 3 for the data {q 0i , q Li }, {s i } and {f ⊥ i , fi } verifying the assumptions (2.7) and (2.11), i = 1, 2, then there exists C > 0 such that The proof of this theorem is based on arguments of monotonicity and will be included in our forthcoming paper [10].Besides the unique weak solvability of Problem 1, it provides the continuous dependence of the solution with respect to the boundary data, the initial gap and the external forces.

Numerical analysis
We now turn to the numerical analysis of the problem.To this end let 0 = L be a partition of the interval [0, L] in N (h) intervals with maximum length h.We denote by Θ h the set of all elements and We consider the finite element spaces where P k (K h ) represents the space of polynomials of degree less or equal than k restricted to K h .It is straightforward to see that V h 3 (0, L) ⊂ H 2 0 (0, L) and 1 (0, L) ⊂ X(0, L), and let P X h : X(0, L) → X h (0, L) denote the projection operator.Then, the discrete version of Problem 3 can be formulated as follows: Using arguments similar to those used in the proof of Theorem 1 it follows that Problem 4 has a unique solution x h = (η h , ϕ h , u h ) ∈ X h (0, L).In addition, the following a priori error estimation holds: where, recall, x = (η, ϕ, u) ∈ X(0, L) is the solution of Problem 3.For a given represents the local Lagrange interpolation operator.Similarly, denote by Π h 3 : C 1 ([0, L]) → V h 3 (0, L) the Hermite global interpolation operator, i.e.
where Π h is the local Hermite interpolation operator.Then, the following interpolation error estimations holds: where v ∈ H r (K), Π h k v denotes its corresponding interpolant in V h k (0, L), l = min{k + 1 − m, r − m} and C > 0 is a constant which does not depend on v, k and h.For the proof of (3.2) see, for instance, [4,6].
Finally, let , where l = min{r 1 , r 2 , r 3 }.Furthermore, by using the previous error estimation and density arguments in (3.1), we conclude that lim h→0 x − x h X(0,L) = 0, which represents a convergence result for our discrete scheme.
Algorithm implementation.The algorithm for the numerical solution of Problem 4 is based on a fixed point strategy to compute the bending, the stretching and the electric potential, iteratively.The novelty lies in the method we use to solve the contact problem on each step.The description of this method represents our main aim in the rest of this section.In order to simplify it, note that the numerical discretization of (2.8) fits in the following general framework: Here and below the symbol E denotes the dual of E. To solve this problem we use a penalization algorithm of the Uzawa family whose performance is improved by combining it with the Newton method.We restrict ourselves to describe the main steps of the algorithm, and refer the reader to [2,11] for further details.Let H ω µ denote the Yosida approximation of the operator H ω = H − ωI, with ωµ < 1.Then the algorithm introduced in [2] is the following: Let q h,0 , ξ h,0 be arbitrary.Given q h,n and ξ h,n−1 , compute q h,n+1 and ξ h,n such that The efficiency of this algorithm is well known, although its convergence is quite slow.For this reason, following [3], we combine it with the Newton method which accelerates its convergence.First, note that if ω = 0, then ωµ < 1 and Therefore, we need to solve the system Besides, given φ ∈ E h , we have where , where Thus, taking p 0 = Bξ h,n+1 + µq h,n+1 and p 1 = Bξ h,n + µq h,n , we obtain that Next, we use the subsets of the mesh Θ h given by and we take into account the values of V in (3.6) to find that which implies that q h,n+1 = 0 on Ω − n and q h,n+1 = k(Bξ h,n+1 −s) on Ω + n .Recall that our aim is to solve the system (3.4)-(3.5).To this end, substituting q h,n+1 in (3.4) by this value obtained by applying the Newton approximation to (3.5), we obtain the following algorithm: Let q h,0 , ξ h,0 be arbitrary.Given q h,n , ξ h,n and Ω + n , compute q h,n+1 , ξ h,n+1 and We finally recall that even though the algorithm above has been presented in a general framework, in our work it is highly simplified since the operator B is, in fact, the identity operator.
In this section we will show the numerical results obtained for several simulations designed to highlight the performance of the algorithm.We consider a beam of length L = 1 m, and a uniformly spaced mesh with element size h = 0.01.The choice of values for the various parameters is the following: Our results are presented in Figures 1-4 where the bending, the stretching, the deformation and the electric potential of the beam are plotted.Next, in Figure 5, we illustrate the influence of the stiffness coefficient on the bending of the beam.We start with the value k = 1 for the stiffness coefficient and increase it up to the value k = 1 × 10 16 .Our results show that more the obstacle is stiff, less the penetration is.We also note that the influence of the obstacle arises only for values of the stiffness coefficient larger than the Young modulus of the beam, E. We also plot the solutions obtained by using both the method in (3.3) and its Newton improvement formulated in (3.8), for k = 1 × 10 16 .As expected, the two solutions are practically the same.Now, in order to show the improvement of the Newton method versus the original algorithm, we represent in Figure 6   to achieve convergence for both algorithms, with various tolerances.As we can see, the convergence of the Newton method is always achieved in about ten iterations, while the Bermúdez-Moreno algorithm (3.3) needs more than one thousand iterations, if the tolerance is smaller than 1 × 10 −3 .In Figures 7, 8 and 9 we show the convergence of the solutions for bending, stretching and electric potential, respectively, as the mesh size decreases.We take the solution obtained for h = 1 × 10 −4 as "exact" solution.We note that the convergence for the electric potential and the stretching field is linear while the convergence for the bending field is slower.
the number of iterations needed