# A Contact Model for Piezoelectric Beams

Á. Rodríguez-Arós, M. T. Cao-Rial, M. Sofonea

► To cite this version:

Á. Rodríguez-Arós, M. T. Cao-Rial, M. Sofonea. A Contact Model for Piezoelectric Beams. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.442-451, 10.1007/978-3-319-55795-3_42. hal-01626911

HAL Id: hal-01626911

https://inria.hal.science/hal-01626911

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A contact model for piezoelectric beams

Á. Rodríguez-Arós¹*, M.T. Cao-Rial², and M. Sofonea³

¹ Departamento de Métodos Matemáticos e Representación, Universidade da Coruña, Paseo de Ronda 51, 15011 A Coruña, Spain

² Technological Institute for Industrial Mathematics (ITMATI), 15782 Campus Vida, Santiago de Compostela, Spain & Departamento de Matemáticas, Universidade da Coruña, Spain

³ Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue de Paul Alduy, 66860 Perpignan, France

**Abstract.** 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.

**Keywords:** Beam, Elasticity, Piezoelectricity, Contact, Normal Compliance, Finite Element Method

## 1 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–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.

## 2 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 $\bar{f}$ and $f^{\perp}$, 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 $\xi : [0, L] \to \mathbb{R}$, a stretching field $u : [0, L] \to \mathbb{R}$ and an electric potential $q : [0, L] \to \mathbb{R}$ such that

$$ (C^{\perp}\xi'')'' = f^{\perp} - p(\xi - s) - \mu_1 R(q) \tilde{R}(\xi - s), \quad (2.1) $$

$$ -(Pu')' - (\varepsilon q')' = \mu_2(\tilde{R}(\xi - s))^2, \quad (2.2) $$

$$ -(\bar{C}u')' + (Pq')' = \bar{f}, \quad (2.3) $$

$$ \xi(0) = \xi'(0) = \xi(L) = \xi'(L) = 0, \quad (2.4) $$

$$ q(0) = q_0, \quad q(L) = q_L, \quad (2.5) $$

$$ u(0) = u(L) = 0. \quad (2.6) $$

Here, $P$ and $\varepsilon > 0$ denote the piezoelectric coefficient and the electric permittivity coefficient, respectively, and $C^{\perp} = EI > 0$, $\bar{C} = EA > 0$. Moreover, $p(\cdot) : \mathbb{R} \to \mathbb{R}_+$ denotes the normal compliance function, which vanishes for negative arguments, and $\mu_1 > 0$, $\mu_2 > 0$ represent coefficients of the system. The dependence of the various functions on $x \in [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 $\tilde{R}(z) = (R(z))_+$, for all $z \in \mathbb{R}$, where $r_+ = \max\{r, 0\}$ denotes the nonnegative part of $r$ and $R$ is the truncation operator given by

$$R(z) = \begin{cases} -M & \text{if } z < -M, \\ z & \text{if } -M \le z \le M, \\ M & \text{if } z > M, \end{cases}$$

$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 $\|\cdot\|_0$, $\|\cdot\|_1$, $\|\cdot\|_2$ and $\|\cdot\|_{\infty}$ the norms on the spaces $L^2(0, L)$, $H^1(0, L)$, $H^2(0, L)$ and $L^{\infty}(0, L)$, respectively. Let $\varphi = q - \hat{q}$, where $\hat{q}$ is a lifting function of $q_0$ and $q_L$ in $H^1(0, L)$ (see [13]). Note that

$$\hat{q} \in C^0([0, L]) \quad \text{and} \quad \max_{x \in [0, L]} \{\hat{q}(x)\} = \max\{q_0, q_L\}.$$

Moreover, without loss of generality, we assume that $q_0, q_L > 0$. As a consequence, it follows that $\hat{q} > 0$ as well. We also assume that

$$s \in C^1([0, L]), \quad s(0) = s'(0) = s(L) = s'(L) = 0, \qquad (2.7)$$

and we denote $\eta = \xi - s$. Then, using a standard procedure, it is easy to obtain the following variational formulation of problem (2.1)–(2.6).

**Problem 2.** Find $\eta \in H_0^2(0, L)$, $u \in H_0^1(0, L)$, $\varphi \in H_0^1(0, L)$, such that

$$\int_0^L C^\perp \eta'' \zeta'' dx + \int_0^L \mu_1 R(\varphi + \hat{q}) \tilde{R}(\eta) \zeta dx + \int_0^L p(\eta) \zeta dx \\ = \int_0^L f^\perp \zeta dx - \int_0^L C^\perp s'' \zeta'' dx, \qquad (2.8)$$

$$\int_0^L P u' \psi' dx + \int_0^L \varepsilon \varphi' \psi' dx - \int_0^L \mu_2 (\tilde{R}(\eta))^2 \psi dx = - \int_0^L \varepsilon \hat{q}' \psi' dx, \qquad (2.9)$$

$$\int_0^L \bar{C} u' v' dx - \int_0^L P \varphi' v' dx = \int_0^L \bar{f} v dx + \int_0^L P \hat{q}' v' dx, \qquad (2.10)$$

for all $\zeta \in H_0^2(0, L)$, $v, \psi \in H_0^1(0, L)$.

In the study of Problem 2 we assume the following hypotheses.

$$\bar{f}, f^\perp \in L^2(0, L), \qquad (2.11)$$

$$C^\perp, P, \varepsilon, \bar{C}, \mu_1, \mu_2 \in L^\infty(0, L), \qquad (2.12)$$

$$C^\perp \ge C_0^\perp > 0, \quad \varepsilon \ge \varepsilon_0 > 0, \quad \bar{C} \ge \bar{C}_0 > 0, \quad \mu_1 > 0, \quad \mu_2 > 0 \text{ a.e. in } (0, L), \qquad (2.13)$$
Moreover, we assume that the normal compliance function $p$ satisfies

$$
\left\{
\begin{array}{ll}
\text{(a)} & p : [0, L] \times \mathbb{R} \to \mathbb{R}_+, \\
\text{(b)} & \text{There exists } c_p > 0 \text{ such that} \\
& |p(x, r_1) - p(x, r_2)| \le c_p |r_1 - r_2|, \\
& \forall r_1, r_2 \in \mathbb{R}, \text{ a.e. } x \in (0, L), \\
\text{(c)} & \text{There exists } m_p \ge 0 \text{ such that} \\
& (p(x, r_1) - p(x, r_2))(r_1 - r_2) \ge m_p |r_1 - r_2|^2, \\
& \forall r_1, r_2 \in \mathbb{R}, \text{ a.e. } x \in (0, L), \\
\text{(d)} & \text{The mapping } p(\cdot, r) : x \mapsto p(x, r) \text{ is measurable on } [0, L], \\
& \text{for all } r \in \mathbb{R}, \\
\text{(e)} & \text{The mapping } p(\cdot, r) : x \mapsto p(x, r) \text{ vanishes for all } r \le 0,
\end{array}
\right.
\qquad (2.14)
$$

and, in addition,

$$
\frac{C_0^{\perp}}{c_2^2 M} > \frac{3}{2} \|\mu_1\|_{\infty} + \|\mu_2\|_{\infty}, \quad \frac{\varepsilon_0}{c_1^2 M} > \frac{1}{2} \|\mu_1\|_{\infty} + \|\mu_2\|_{\infty}. \qquad (2.15)
$$

Next, we consider the space $X(0, L) = H_0^2(0, L) \times H_0^1(0, L) \times H_0^1(0, L)$, which is a real Hilbert space with the canonical inner product, denoted by $(\cdot, \cdot)_{X(0,L)}$. We then define the operator $A : X(0, L) \to X(0, L)$ and the element $\mathbf{F} \in X(0, L)$ by equalities

$$
\begin{equation}
\begin{split}
(A\mathbf{x}, \mathbf{y})_{X(0,L)} = \int_0^L C^\perp \eta'' \zeta'' dx &+ \int_0^L \mu_1 R(\varphi + \hat{q}) \tilde{R}(\eta) \zeta dx + \int_0^L p(\eta) \zeta dx \\
&+ \int_0^L P u' \psi' dx + \int_0^L \varepsilon \varphi' \psi' dx - \int_0^L \mu_2 (\tilde{R}(\eta))^2 \psi dx + \int_0^L \bar{C} u' v' dx \\
&- \int_0^L P \varphi' v' dx \quad \forall \mathbf{x} = (\eta, \varphi, u), \mathbf{y} = (\zeta, \psi, v) \in X(0, L),
\end{split}
\tag{2.16}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
(F, y)_{X(0,L)} = \int_0^L f^{\perp} \zeta dx &+ \int_0^L \bar{f} v dx - \int_0^L C^{\perp} s'' \zeta'' dx - \int_0^L \varepsilon \hat{q}' \psi' dx \\
&+ \int_0^L P \hat{q}' v' dx \quad \forall \mathbf{y} = (\zeta, \psi, v) \in X(0, L).
\end{split}
\tag{2.17}
\end{equation}
$$

Then, an equivalent formulation of Problem 2 is as follows.

Problem 3. Find $\mathbf{x} = (\eta, \varphi, u) \in X(0, L)$ such that

$$
(A\mathbf{x}, \mathbf{y}) = (\mathbf{F}, \mathbf{y})_{X(0,L)} \quad \forall \mathbf{y} \in X(0, L). \qquad (2.18)
$$

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 $\mathbf{x} = (\eta, \varphi, u) \in X(0, L)$ to Problem 3. Moreover, if $\mathbf{x}_i = (\eta_i, u_i, \varphi_i)$
represents the solution of Problem 3 for the data {$q_{0i}, q_{Li}$}, {$s_i$} and {$f_i^\perp, \bar{f}_i$} verifying the assumptions (2.7) and (2.11), $i = 1, 2$, then there exists $C > 0$ such that

$$ \|\mathbf{x}_1 - \mathbf{x}_2\|_{X(0,L)} \le C(\|f_1^\perp - f_2^\perp\|_0 + \|s_1 - s_2\|_2 + \|\bar{f}_1 - \bar{f}_2\|_0 + \|\hat{q}_1 - \hat{q}_2\|_1). \quad (2.19) $$

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.

### 3 Numerical analysis

We now turn to the numerical analysis of the problem. To this end let $0 < h < L$, $N(h) \in \mathbb{N}$ and let $0 = x_0^h < x_1^h < \dots < x_i^h < x_{i+1}^h < \dots < x_{N(h)}^h = L$ be a partition of the interval $[0, L]$ in $N(h)$ intervals with maximum length $h$. We denote by $\Theta^h$ the set of all elements and $K_i^h = [x_i^h, x_{i+1}^h] \in \Theta^h$, $0 \le i < N(h)$. We consider the finite element spaces

$$ \begin{aligned} V_1^h(0, L) &= \{\xi^h \in C^0([0, L]), \xi_{|K_i^h}^h \in P_1(K_i^h), 0 \le i < N(h), \xi^h(0) = \xi^h(L) = 0\}, \\ V_3^h(0, L) &= \{\xi^h \in C^1([0, L]), \xi_{|K_i^h}^h \in P_3(K_i^h), 0 \le i < N(h), \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \xi^h(0) = \xi^h(L) = (\xi^h)'(0) = (\xi^h)'(L) = 0\}. \end{aligned} $$

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_3^h(0, L) \subset H_0^2(0, L)$ and $V_1^h(0, L) \subset H_0^1(0, L)$. Let $X^h(0, L) = V_3^h(0, L) \times V_1^h(0, L) \times V_1^h(0, L) \subset X(0, L)$, and let $P_{X^h}: X(0, L) \to X^h(0, L)$ denote the projection operator. Then, the discrete version of Problem 3 can be formulated as follows:

**Problem 4.** Find $\mathbf{x}^h = (\eta^h, \varphi^h, u^h) \in X^h(0, L)$ such that

$$ (A\mathbf{x}^h, \mathbf{y}^h)_{X(0,L)} = (\mathbf{F}, \mathbf{y}^h)_{X(0,L)} \quad \forall \mathbf{y}^h = (\zeta^h, \psi^h, v^h) \in X^h(0,L). $$

Using arguments similar to those used in the proof of Theorem 1 it follows that Problem 4 has a unique solution $\mathbf{x}^h = (\eta^h, \varphi^h, u^h) \in X^h(0, L)$. In addition, the following a priori error estimation holds:

$$ \|\mathbf{x} - \mathbf{x}^h\|_{X(0,L)} \le C \|\mathbf{x} - \mathbf{y}^h\|_{X(0,L)} \quad \forall \mathbf{y}^h \in X^h(0,L), \quad (3.1) $$

where, recall, $\mathbf{x} = (\eta, \varphi, u) \in X(0, L)$ is the solution of Problem 3.

For a given $K_i^h = [x_i^h, x_{i+1}^h] \in \Theta^h$, we use a local nodal notation, so $K_i^h = [x_0^{i,h}, x_1^{i,h}]$. Denote by $\Pi_1^h: C([0, L]) \to V_1^h(0, L)$ the Lagrange global interpolation operator, i.e.

$$ \Pi_1^h v|_{K_i^h} = \Pi_1^h v|_{K_i^h} v \quad \forall K_i^h \in \Theta^h, v \in C([0, L]), $$
where $\Pi_{1|K_i^h}^h : C(K_i^h) \to P_1(K_i^h)$ represents the local Lagrange interpolation operator. Similarly, denote by $\Pi_3^h : C^1([0, L]) \to V_3^h(0, L)$ the Hermite global interpolation operator, i.e.

$$ \Pi_3^h v|_{K_i^h} = \Pi_{3|K_i^h}^h v \quad \forall K_i^h \in \Theta^h \quad v \in C^1([0, L]), $$

where $\Pi_{3|K_i^h}^h : C^1(K_i^h) \to P_3(K_i^h)$ is the local Hermite interpolation operator. Then, the following interpolation error estimations holds:

$$ \|v - \Pi_k^h v\|_{m,K} \le C h_K^l |v|_{r,K} \quad \forall K \in \Theta^h, \quad (3.2) $$

where $v \in H^r(K)$, $\Pi_k^h v$ denotes its corresponding interpolant in $V_k^h(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 $\Pi^h : C^1([0, L]) \times C^0([0, L]) \times C^0([0, L]) \to X^h(0, L)$ denote the global interpolation operator given by $\Pi^h(y) = (\Pi_3^h(\zeta), \Pi_1^h(\psi), \Pi_1^h(v))$ for all $y = (\zeta, \psi, v) \in C^1([0, L]) \times C^0([0, L]) \times C^0([0, L])$. Therefore, given $\{r_1, r_2, r_3\} \subset \mathbb{N}$, from (3.2) we find that

$$ \|y - \Pi^h y\|_{X(0,L)} \le C (h^{r_1} \|\zeta\|_2 + h^{r_2} \|\psi\|_1 + h^{r_3} \|v\|_1) \le C h^l \|y\|_{X(0,L)}, $$

for all $y = (\zeta, \psi, v) \in H^{2+r_1}(0, L) \times H^{1+r_2}(0, L) \times H^{1+r_3}(0, L)$, 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 \to 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:

Find $\xi^h \in E^h$ such that $L \in A\xi^h + BHB^*\xi^h$, where $H : E^h \to 2^{E^h}$ is a maximal monotone operator, $B \in \mathcal{L}(E^h, V_3^{h'})$, $A \in \mathcal{L}(V_3^{h}, V_3^{h'})$ and $L \in V_3^{h'}$.

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_\mu^\omega$ denote the Yosida approximation of the operator $H^\omega = H - \omega I$, with $\omega\mu < 1$. Then the algorithm introduced in [2] is the following:

Let $q^{h,0}, \xi^{h,0}$ be arbitrary. Given $q^{h,n}$ and $\xi^{h,n-1}$, compute $q^{h,n+1}$ and $\xi^{h,n}$ such that

$$ \left\{ \begin{array}{l} A\xi^{h,n} + \omega BB^*\xi^{h,n} = L - Bq^{h,n}, \\ q^{h,n+1/2} = H_\mu^\omega [B^*\xi^{h,n} + \mu q^{h,n}], \\ q^{h,n+1} = \rho q^{h,n+1/2} + (1-\rho)q^{h,n}. \end{array} \right. \qquad (3.3) $$
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 $\omega = 0$, then $\omega\mu < 1$ and

$$H_{\mu}^{0}(x) = H_{\mu}(x) = \frac{k}{1 + k\mu}(I - \Pi_{E^{h}})(x).$$

Therefore, we need to solve the system

$$A\xi^{h,n+1} = L - Bq^{h,n+1}, \quad (3.4)$$

$$q^{h,n+1} = H_{\mu}(B\xi^{h,n+1} + \mu q^{h,n+1}). \quad (3.5)$$

Besides, given $\phi \in E^h$, we have $H_\mu(\phi)|_C = G(\phi|_C)$ where

$$G(p) = \begin{cases} 0 & \text{if } p \le s \\ \frac{k}{1+k\mu}(p-s)(x) & \text{if } p > s. \end{cases}$$

We approximate $G(p_0)$ by $G(p_1) + V(p_0 - p_1) + O(|p_0 - p_1|^2)$, where

$$V \in \partial G(p_0) = \begin{cases} 0 & \text{if } p_0 < s \\ \left[ 0, \frac{k}{1+k\mu} \right] & \text{if } p_0 = s. \\ \frac{k}{1+k\mu} & \text{if } p_0 > s \end{cases}. \quad (3.6)$$

Thus, taking $p_0 = B\xi^{h,n+1} + \mu q^{h,n+1}$ and $p_1 = B\xi^{h,n} + \mu q^{h,n}$, we obtain that
$$q^{h,n+1} = H_{\mu}(p_0) = G(B\xi^{h,n} + \mu q^{h,n}) + V(B\xi^{h,n+1} + \mu q^{h,n+1} - B\xi^{h,n} + \mu q^{h,n}).$$

Next, we use the subsets of the mesh $\Theta^h$ given by

$$\Omega_n^+ = \{ K_i^h \in \Theta^h ; (B\xi^{h,n} + \mu q^{h,n})|_{K_i^h} > s \}, \quad \Omega_n^- = \Theta^h \setminus \Omega_n^+; \quad (3.7)$$

and we take into account the values of $V$ in (3.6) to find that

$$q^{h,n+1} = \begin{cases} 0 & \text{on } \Omega_n^-, \\ \frac{k}{1+k\mu}(B\xi^{h,n+1} + \mu q^{h,n+1}) - \frac{k}{1+k\mu}s & \text{on } \Omega_n^+, \end{cases}$$

which implies that $q^{h,n+1} = 0$ on $\Omega_n^-$ and $q^{h,n+1} = k(B\xi^{h,n+1} - s)$ on $\Omega_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}$, $\xi^{h,0}$ be arbitrary. Given $q^{h,n}$, $\xi^{h,n}$ and $\Omega_n^+$, compute $q^{h,n+1}$, $\xi^{h,n+1}$ and $\Omega_{n+1}^+$ such that

$$\left\{ \begin{array}{l} \displaystyle \int_{\Omega} A \xi^{h,n+1} \zeta dx + \int_{\Omega_n^+} k B \xi^{h,n+1} \zeta dx = \int_{\Omega} L \zeta dx + \int_{\Omega_n^+} k s \zeta dx, \\[1.5ex] \displaystyle q^{h,n+1} = \begin{cases} 0 & \text{on } \Omega_n^-, \\ k(B\xi^{h,n+1} - s) & \text{on } \Omega_n^+, \end{cases} \\[1.5ex] \displaystyle \Omega_{n+1}^+ = \{ K_i^h \in \Theta^h ; (B\xi^{h,n+1} + \mu q^{h,n+1})|_{K_i^h} > s \}. \end{array} \right. \quad (3.8)$$

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.
## 4 Numerical simulations

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:

$$E = 1 \times 10^5 \frac{N}{m^2}, I = 0.05 \, m^4, A = 1 \, m^2, P = 100 \frac{N \cdot m}{V},$$

$$\mu_1 = \mu_2 = 1 \frac{N}{V \cdot m^2}, \varepsilon = 1 \frac{N \cdot m^2}{V}, \bar{f} = 1 \times 10^6 \frac{N}{n}, f^{\perp} = 1 \times 10^6 \frac{N}{m}.$$

Our results are presented in Figures 1–4 where the bending, the stretching, the deformation and the electric potential of the beam are plotted.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Bending.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Stretching.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Deformation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Electric potential.</figcaption></figure>

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 \times 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 \times 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 the number of iterations needed
<figure><img src="image_6.png" /><figcaption>Fig. 5. Bending for several coefficients of stiffness</figcaption></figure>

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 \times 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 \times 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.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Iterations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Bending Error.</figcaption></figure>

References

1. Bellis, S., Imperiale, S.: Dynamical 1D models of passive piezoelectric sensors. Math. Mech. Solids (2013).

2. Bermúdez, A., Moreno, C.: Duality methods for solving variational inequalities. Comput. Math. Appl. 7(1), 43–58 (1981)

3. Cao, M.T., Moreno, C., Quintela, P.: Simulation of Rayleigh waves in cracked plates. Mathematical Methods in the Applied Sciences 30, 15–42 (2007).
<figure><img src="image_9.png" /><figcaption>Fig. 8. Stretching Error</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Electric Potential Error</figcaption></figure>

4. Ciarlet, P.G.: The finite element method for elliptic problems. Studies in Mathematics and its Applications, Vol. 4. North-Holland Publishing Co., Amsterdam-New York-Oxford (1978)

5. Figueiredo, I., Leal, C.: A generalized piezoelectric Bernoulli-Navier anisotropic rod model. J. Elasticity 31, 85–106 (2006).

6. Hughes, T.J.R.: The finite element method. Prentice Hall, Inc., Englewood Cliffs, NJ (1987)

7. Irago, H., Viaño, J.M., Rodríguez-Arós, Á.: Asymptotic derivation of frictionless contact models for elastic rods on a foundation with normal compliance. Nonlinear Anal. Real World Appl. 14, 852–866 (2013).

8. Ribeiro, C.: Asymptotic derivation of models for anisotropic piezoelectric beams and shallow arches. Ph.D. Thesis, Department of Mathematics and Applications, Minho University (2009).

9. Rodríguez-Arós, Á., Viaño, J.M.: A bending-stretching model in adhesive contact for elastic rods obtained by using asymptotic methods. Nonlinear Analysis-RWA 22, 632–644 (2015).

10. Rodríguez-Arós, Á., Cao-Rial, M.T., Sofonea, M.: Mathematical analysis of a piezoelectric elastic beam in contact with a deformable foundation. Preprint.

11. Viaño, J.M.: Análisis numérico de un método numérico con elementos finitos para problemas de contacto unilateral sin rozamiento en lasticidad: Aproximación y resolución de los problemas discretos. Rev. Internac. Métod. Numér. Cálc. Diseñ. Ingr. 2, 63–86 (1986)

12. Viaño, J.M.: The one-dimensional obstacle problem as approximation of the three-dimensional Signorini problem. Bull. Math. Soc. Sci. Math. Roumanie (N.S.) 48, 243–258 (2005).

13. Viaño, J., Figueiredo, J., Ribeiro, C., Rodríguez-Arós, Á.: A model for bending and stretching of piezoelectric rods obtained by asymptotic analysis. Zeitschrift fr angewandte Mathematik und Physik 66(3), 1207–1232 (2015)

14. Viaño, J.M., Ribeiro, C., Figueiredo, J.: Asymptotic modelling of a piezoelectric beam. In: Proceedings of Congreso de Métodos Numéricos en Ingeniería, pp. 1–17. Granada, Spain, A310, (2005).

15. Viaño, J.M., Ribeiro, C., Figueiredo, J.: Asymptotic modelling of a piezoelectric beam. In: Proceedings of II ECCOMAS Thematic Conference of Smart Materials and Structures, pp. 1–12. Lisbon, Portugal, EO26MOD (2005),

16. Viaño, J.M., Rodríguez-Arós, Á., Sofonea, M.: Asymptotic derivation of quasistatic frictional contact models with wear for elastic rods. J. Math. Anal. Appl. 401, 641–653 (2013).

17. Weller, T., Licht, C.: Asymptotic modeling of linearly piezoelectric slender rods. C. R. Mecanique 336, 572–577 (2008).