3D BEM for orthotropic frictional contact of piezoelectric bodies

Abstract

A numerical methodology to model the three-dimensional frictional contact interaction of piezoelectric materials in presence of electric fields is presented in this work. The boundary element method (BEM) is used in order to compute the electro-elastic influence coefficients. The proposed BEM formulation employs an explicit approach for the evaluation of the corresponding fundamental solutions, which are valid for general anisotropic behaviour meanwhile mathematical degeneracies in the context of the Stroh formalism are allowed. The contact methodology is based on an augmented Lagrangian formulation and uses an iterative Uzawa scheme of resolution. An orthotropic frictional law is implemented in this work so anisotropy is present both in the bulk and in the surface. The methodology is validated by comparison with benchmark analytical solutions. Some additional examples are presented and discussed in detail, revealing the importance of considering orthotropic frictional contact conditions in the electro-elastic analysis of this kind of problems.

Appendices

Appendix 1: Barnett & Lothe representation

As proposed by Barnett & Lothe [53], the linear PE problem may be formulated in an elastic-like fashion by considering a generalized displacement vector extended with the electric potential as

$$\begin{aligned} u_{J}=\left\{ \begin{array} [c]{ll}u_{j} &{} \quad J\leqslant 3\\ \varphi &{} \quad J=4, \end{array} \right. \end{aligned}$$

(35)

a traction vector extended with electric charge flux as

$$\begin{aligned} t_{J}=\left\{ \begin{array} [c]{ll} t_{j} &{} \quad J\leqslant 3\\ q &{} \quad J=4, \end{array} \right. \end{aligned}$$

(36)

a stress tensor extended with the electric displacements as

$$\begin{aligned} \sigma _{iJ}=\left\{ \begin{array} [c]{ll} \sigma _{ij} &{} \quad J\leqslant 3\\ D_{i} &{} \quad J=4, \end{array} \right. \end{aligned}$$

(37)

and an extended elasticity tensor with the following components

$$\begin{aligned} C_{iJKm}=\left\{ \begin{array} [c]{ll} c_{ijkm} &{} \quad J,K\leqslant 3\\ e_{mij} &{} \quad J\leqslant 3;\;K=4\\ e_{ikm} &{} \quad J=4;\;K\leqslant 3\\ -\epsilon _{im} &{} \quad J,K=4. \end{array} \right. \end{aligned}$$

(38)

Symmetry $C_{iJKm}=C_{mKJi}$ is satisfied by virtue of (4). In the above definitions the lowercase (elastic) and uppercase (extended) subscripts take values 1, 2, 3 and 1, 2, 3 (elastic), 4 (electric), respectively. Then, with this compact representation the extended equilibrium Eq. (1) can be expressed as

$$\begin{aligned} \sigma _{iJ,i}=0 \end{aligned}$$

(39)

and the constitutive Eqs. (3) as

$$\begin{aligned} \sigma _{iJ}=C_{iJKm}u_{K,m}. \end{aligned}$$

(40)

So the Barnett & Lothe representation allows also to arrive to a PE BEM formulation in a very compact way.

Appendix 2: Fundamental solutions

The explicit fundamental solutions recently proposed by Buroni & Sáez [36] have been implemented in the presented formulation and they are briefly described in this appendix. In homogeneous media they depend on the relative vector $\mathbf {x}-\mathbf {x'}$, henceforth, for simplicity, it is considered that the Cartesian coordinate system $(x_i)$ ($i=1,2,3$) has the origin at the collocation point $\mathbf {x'}$.

Extended displacement fundamental solution can be expressed as a singular term by a modulation function $\mathbf {H}$ as

$$\begin{aligned} \check{U}_{JK}(\mathbf x )=\frac{1}{4\pi r}H_{JK}(\mathbf x ) \end{aligned}$$

(41)

where $\mathbf {x}=r\hat{\mathbf {e}}$ with $r=|\mathbf {x}|\ne 0$. The modulation function $H_{JK}(\mathbf {x})$ depends on the direction of $\mathbf {x}$ but not on its modulus, so $H_{JK}(\mathbf {x})=H_{JK}(\hat{\mathbf {e}})$ and that is one of the three extended Barnett-Lothe tensors which is symmetric and $\mathbf {H}(\hat{\mathbf {e}})=\mathbf {H}(\hat{\mathbf {-e}})$. Hence, $\check{\mathbf {U}}(\mathbf {x})$ is also symmetric and $\check{\mathbf {U}}(\mathbf {x})=\check{\mathbf {U}}(\mathbf {-x})$. The tensor $H_{JK}$ can be evaluated as

$$\begin{aligned} H_{JK}(\hat{\mathbf {e}})=\frac{1}{\pi }\int _{-\infty }^{+\infty }\varGamma _{JK}^{-1}(p)dp, \end{aligned}$$

(42)

with

$$\begin{aligned} \varGamma _{JK}(p)=Q_{JK}+(R_{JK}+R_{KJ})p+T_{JK}p^2, \end{aligned}$$

(43)

being

$$\begin{aligned}&Q_{JK}=C_{iJKm}n_in_m,\quad R_{JK}=C_{iJKm}n_im_m,\quad \nonumber \\&T_{JK}=C_{iJKm}m_im_m, \end{aligned}$$

(44)

where $n_i$ and $m_i$ are the components of any two mutually orthogonal unit vectors such that $(\mathbf {n},\mathbf {m},\hat{\mathbf {e}})$ is a right-handed triad. Note that $Q_{JK}$ and $T_{JK}$ are symmetric like their elastic counterparts [54], but the PE coupling induces a loss of positive definiteness of these matrices. However, $Q_{JK}$, $T_{JK}$ and $\varGamma _{JK}$ are non-singular so that their inverses are guaranteed. Moreover, $\mathbf {H}$ and $\check{\mathbf {U}}$ are independent of the choice of the unit vectors $\mathbf {m}$ and $\mathbf {n}$ on the oblique plane. The kernel in Eq. (42) is a single-valued holomorphic function in the upper complex half-plane except at four complex poles with positive imaginary part and their conjugates that corresponds to the roots of the eight-order polynomial characteristic equation

$$\begin{aligned} |\varvec{\varGamma }(p)|=0. \end{aligned}$$

(45)

The determinant in (45) can be factorized as

$$\begin{aligned} |\varvec{\varGamma }(p)|=|\mathbf {T}| \prod _{\xi =1}^4 (p-p_\xi )(p-\bar{p}_\xi ), \end{aligned}$$

(46)

where $p_\xi$ are known as the Stroh’s eigenvalues and the bar over $p_\xi$ denotes complex conjugate and $\mathbf {T}$ as defined in Eq. (44). Therefore, the integration in Eq.  (42) can be done by the Cauchy’s residue theory. At most, there are N ($1\leqq N\leqq 4$) distinct Stroh’s eigenvalues $p_\alpha$ of $m _\alpha$-multiplicity. Hence, a general expression, valid for degenerate and non-degenerate materials, of the extended Barnett-Lothe tensor $H_{JK}$ is obtained as [36]

$$\begin{aligned}&H_{JK}(\hat{\mathbf {e}})=\frac{2 \mathrm {i}}{|\mathbf {T}|}\sum _{\alpha =1}^N\frac{1}{(m _\alpha -1)!}\nonumber \\&\quad \left[ \frac{d^{m _\alpha -1}}{dp^{m _\alpha -1}} \left\{ \!\frac{\hat{\varGamma }_{JK}(p)}{(p-\bar{p}_\alpha )^{m_\alpha } \prod \limits _{\begin{array}{c} \xi =1\\ \xi \ne \alpha \end{array}}^N [(p-p_\xi )(p-\bar{p}_\xi )]^{m_\xi }} \right\} \!\right] _{at\,p=p_\alpha }.\nonumber \\ \end{aligned}$$

(47)

where $\hat{\varGamma }_{JK}$ is the adjoint of $\varGamma _{JK}$ defined as $\varGamma _{PJ}(p)\hat{\varGamma }_{JK}(p)\!=|\varvec{\varGamma }(p)|\delta _{PK}$ and $\mathrm {i}=\sqrt{-1}$.

The extended traction fundamental solution follows from the derivative of the extended displacement fundamental solution as

$$\begin{aligned} \check{T}_{JK}=C_{iJMl}\check{U}_{MK,l}\nu _i \end{aligned}$$

(48)

where $\nu _i$ are the components of the external unit normal vector to the boundary $\partial \varOmega$ at point $\mathbf {x}$. In similar way to Eq. (41), the derivative of the displacement fundamental solution may be expressed as

$$\begin{aligned} \check{U}_{PJ,q}(\mathbf {x})=\frac{1}{4 \pi r^2}\tilde{U}_{PJq}(\hat{\mathbf {e}}) \end{aligned}$$

(49)

where the modulation function is

$$\begin{aligned} \tilde{U}_{PJq}(\hat{\mathbf {e}})= & {} -\hat{e}_{q}H_{PJ}\nonumber \\&+\frac{C_{rKMs}}{\pi }\left( M_{qsPKMJ}\hat{e}_r+M_{qrPKMJ}\hat{e}_s\right) , \end{aligned}$$

(50)

that only depends on the orientation of $\mathbf {x}$ ($\hat{\mathbf {e}}$) but not on its modulus r. The $M_{ijPKMN}$ components have the following integral representation in terms of the parameter p

$$\begin{aligned}&M_{ijPKMN}(\hat{\mathbf {e}})\nonumber \\&\quad =\frac{1}{|\mathbf {T}|^2}\int _{-\infty }^{+\infty }\frac{\varPhi _{ijPKMN}(p)}{(p-{p}_1)^2(p-{p}_2)^2(p-{p}_3)^2(p-{p}_4)^2}dp,\nonumber \\ \end{aligned}$$

(51)

where $\mathbf {T}$ has been previously defined in (44), $p_\alpha$ are the Stroh’s eigenvalues and the function

$$\begin{aligned} \varPhi _{ijPKMN}(p):=\frac{\tilde{B}_{ij}(p)\hat{\varGamma }_{PK}(p)\hat{\varGamma }_{MN}(p)}{(p-\bar{p}_1)^2(p-\bar{p}_2)^2(p-\bar{p}_3)^2(p-\bar{p}_4)^2} \end{aligned}$$

(52)

has been introduced together with definition

$$\begin{aligned} \tilde{B}_{ij}:=n_in_j+(n_im_j+m_in_j)p+m_im_jp^2. \end{aligned}$$

(53)

$\varPhi _{ijPKMN}(p)$ is a holomorphic function everywhere in the upper half-plane and the kernel in the integral (51) has four complex-double poles with positive imaginary part corresponding to the roots of $|\varvec{\varGamma }(p)|^2=0$ (no poles on the real axis). Hence, like in the integration of Eq. (42), a general expression both for degenerate and non-degenerate materials is derived for the $M_{ijPKMN}$ components to yield [36]

$$\begin{aligned}&M_{ijPKMN}(\hat{\mathbf {e}})=\frac{2\pi \mathrm {i}}{|\mathbf {T}|^2}\sum _{\alpha =1}^N\frac{1}{(2m _\alpha -1)!}\nonumber \\&\quad \left[ \frac{d^{2m _\alpha -1}}{dp^{2m _\alpha -1}}\left\{ \frac{\varPhi _{ijPKMN}(p)}{\prod \limits _{\begin{array}{c} \xi =1\\ \xi \ne \alpha \end{array}}^N [(p-p_\xi )]^{2m_\xi }}\right\} \right] _{at\,p=p_\alpha }. \end{aligned}$$

(54)

Because of the symmetry of $\tilde{B}_{ij}$ and the adjoint matrix $\varvec{\hat{\varGamma }}$, the components $M_{qsPKMJ}$ satisfy the following symmetry conditions

$$\begin{aligned} M_{qsPKMJ}= & {} M_{qsMJPK}=M_{qsKPMJ}=M_{qsPKJM}, \end{aligned}$$

(55)

$$\begin{aligned} M_{qsPKMJ}= & {} M_{sqPKMJ}. \end{aligned}$$

(56)

These symmetries allow considerable reduction in the number of components $M_{qsPKMJ}$ to be calculated, and must be considered in the numerical implementation.