# A regularized direct symmetric Galerkin BIE formulation for three-dimensional elastoplasticity

Bertrand Burgardt, Marc Bonnet, Anh Le Van

► To cite this version:

Bertrand Burgardt, Marc Bonnet, Anh Le Van. A regularized direct symmetric Galerkin BIE formulation for three-dimensional elastoplasticity. Symposium on Advanced Mathematical and Computational Mechanics Aspects of the Boundary Element Method, May 1999, Cracovie, Poland. pp.51-61, 10.1007/978-94-015-9793-7_5 . hal-01002840

HAL Id: hal-01002840

https://hal.science/hal-01002840

Submitted on 13 Oct 2016

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 REGULARIZED DIRECT SYMMETRIC GALERKIN BIE FORMULATION FOR THREE-DIMENSIONAL ELASTOPLASTICITY

Bertrand Burgardt(a), Marc Bonnet(b), Anh Le Van(a)

Bertrand.Burgardt@ec-nantes.fr, bonnet@lms.polytechnique.fr, Anh.Le-Van@ec-nantes.fr

**Abstract** This paper deals with a symmetric regularized variational boundary/domain formulation for quasi-static 3D elastoplasticity, which is shown to express the stationarity of a certain energy functional. An implicit constitutive integration scheme is implemented. The global consistent tangent operator associated with the non-linear algebraic system of equations is shown to be symmetric.

**Keywords:** Symmetric Galerkin BEM, elastoplasticity, implicit constitutive integration

## 1. THE QUASI-STATIC ELASTOPLASTIC MODEL

Mechanical problems with material or geometrical non-linearities can be solved with BEM. In that case, the boundary integrals are supplemented with domain integrals containing plastic strains treated as unknown distributions of initial strains. Non-linear domain-BEM formulations are attractive in some specific situations like infinite media, fracture mechanics problems, ... This paper addresses the formulation and numerical implementation of a symmetric Galerkin boundary integral equation (SGBIE) method for solving three-dimensional small-strain quasi-static elastic-plastic problems.

Consider a homogeneous body occupying the open domain Ω, referred to a Cartesian orthogonal system. Its piecewise smooth boundary S is split into two disjoint open subsets $S_u$ and $S_t$, over which histories of boundary data (displacement $u^0(\cdot, t)$ and traction $t^0(\cdot, t)$, respectively) are prescribed in a quasi-static manner, $t$ denoting a monotonically in-

<footnote>(a) Laboratoire Mécanique et Matériaux, Ecole Centrale de Nantes, France</footnote>

<footnote>(b) Laboratoire de Mecanique des Solides (CNRS UMR 7649), Ecole Polytechnique, France</footnote>
creasing time-like parameter hereafter referred to simply as ‘time’. Let $\Omega_p$ denote the potentially plastic region, i.e. an open bounded subset of $\Omega$ outside of which plastic stains $\varepsilon^p$ and internal variables $\alpha$ are assumed to vanish.

Assuming small displacements and strains, the quasi-static evolution of the considered solid from the initial state is governed (see [9]) by (i) mechanical field equations (equilibrium, compatibility) and boundary conditions, (ii) state laws (generalized Hooke's laws) and (iii) evolution laws (plastic flow rules, consistency condition). The latter equations imply that any numerical solution procedure should involve a time-stepping scheme. In this paper, a single-step integration method [21], namely the well-known backward difference method, is used; it consists in solving an increment problem for every time step with enforcement of the constitutive equations at the step end. This procedure is nowadays quite popular and some interesting results about consistency, convergence and stability are known [18]. After discretization in time, the equations of the evolutive problem are:

$$ (\text{field equations}) \qquad \operatorname{div} \boldsymbol{\sigma}_{n+1} = \mathbf{0} \tag{1} $$

$$ (\text{boundary cond.}) \qquad \begin{aligned} \boldsymbol{u}_{n+1} &= \boldsymbol{u}_{n+1}^0 & \text{on } S_u \\ \boldsymbol{\sigma}_{n+1} \cdot \boldsymbol{n} &= t_{n+1}^0 & \text{on } S_t \end{aligned} \tag{2} $$

$$ (\text{state laws}) \qquad \begin{aligned} \boldsymbol{\sigma}_{n+1} &= C : (\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_{n+1}^p) \\ \boldsymbol{q}_{n+1} &= \frac{\partial \Theta}{\partial \boldsymbol{\alpha}} \end{aligned} \tag{3} $$

$$ (\text{evolution laws}) \qquad \begin{aligned} \boldsymbol{\varepsilon}_{n+1}^p &= \boldsymbol{\varepsilon}_n^p + \lambda_{n+1} \frac{\partial \phi}{\partial \boldsymbol{\sigma}}(\boldsymbol{\sigma}_{n+1}, \boldsymbol{q}_{n+1}) \\ \boldsymbol{\alpha}_{n+1} &= \boldsymbol{\alpha}_n - \lambda_{n+1} \frac{\partial \phi}{\partial \boldsymbol{q}}(\boldsymbol{\sigma}_{n+1}, \boldsymbol{q}_{n+1}) \end{aligned} \tag{4} $$
with $\lambda_{n+1} \phi(\boldsymbol{\sigma}_{n+1}, \boldsymbol{q}_{n+1}) = 0$
$\lambda_{n+1} \ge 0 \quad \phi(\boldsymbol{\sigma}_{n+1}, \boldsymbol{q}_{n+1}) \le 0$

where $\boldsymbol{\sigma}$ and $\varepsilon$ denote stresses and total strains, $\boldsymbol{n}$ the unit normal to $S$ outward to $\Omega$, $C$ the elastic stiffness tensor, $\boldsymbol{q}$ the thermodynamical force associated to $\boldsymbol{\alpha}$, $\Theta$ the hardening potential and $\phi$ the yield surface.

This set of relations can be solved by a FEM scheme (based on variational principles in either displacement [20], mixed [19] or Hu-Washizu [21] forms), but also by means of boundary/domain formulations.
## 2. VARIATIONAL INTEGRAL FORMULATION

Boundary element methods for quasi-static elastoplasticity based on variational principles were proposed in the late eighties [12, 14] and have since seen several implementations for two-dimensional problems [3]. Weakly singular SGBIE formulations for three-dimensional elasticity are expounded in [2] (as stationarity equations of the elastic potential energy for displacement fields in elastic equilibrium) and [10]. The former formulation is here extended to small-strain elastoplasticity.

Let $\zeta_{n+1} = (u_{n+1}, \sigma_{n+1}, \varepsilon_{n+1}, \varepsilon_{n+1}^p, q_{n+1}, \alpha_{n+1})$ denote the solution to the governing relations at time $t_{n+1}$. The following representation formula holds for any displacement solving the field equations (1) and Hooke's law:

$$ u_k(\boldsymbol{x}) = \int_S \tilde{t}_i(\boldsymbol{y}) U_i^k(\boldsymbol{y}-\boldsymbol{x}) \, dS_y - \int_S \sum_{ij}^k (\boldsymbol{y}-\boldsymbol{x}) n_j(\boldsymbol{y}) \tilde{u}_i(\boldsymbol{y}) \, dS_y + \int_{\Omega} \sum_{ij}^k (\boldsymbol{y}-\boldsymbol{x}) \varepsilon_{ij}^p(\boldsymbol{y}) \, dV_y \quad (5) $$

where $\tilde{u}, \tilde{t}$ are the boundary displacement and traction associated with $u$, respectively, $U(r)$ and $\Sigma(r) = C : \nabla U(r)$ the Kelvin fundamental displacement and elastic stress, respectively; strains and stresses at $x$ are then obtained upon differentiation of the above formula with respect to $x$, considering $\varepsilon^p$ as initial strains.

Following [2] and [16], the variational equation on the submanifold $\chi_{n+1}^\sigma = ((\tilde{u}, \tilde{t}), \sigma, \varepsilon^p, q, \alpha)_{n+1}$ of $\zeta_{n+1}$.

$$ \begin{aligned} & \int_{S_t} \tilde{u} \cdot [C : \nabla^s \delta \boldsymbol{u}] \cdot \boldsymbol{n} \, dS + \int_{S_u} \tilde{u}^0 \cdot [C : \nabla^s \delta \boldsymbol{u}] \cdot \boldsymbol{n} \, dS \\ & - \int_{S_t} \tilde{t}^0 \cdot \gamma_0 \delta \boldsymbol{u} \, dS - \int_{S_u} \tilde{t} \cdot \gamma_0 \delta \boldsymbol{u} \, dS - \int_{\Omega_p} \nabla \delta \boldsymbol{u} : C : \varepsilon^p \, dV \\ & + \int_{\Omega_p} \delta \alpha \cdot \left[ \frac{\partial \Theta}{\partial \alpha}(\boldsymbol{\alpha}) - \boldsymbol{q} \right] \, dV + \int_{\Omega_p} \delta \varepsilon^p : [\boldsymbol{\sigma} - C : (\nabla^s \boldsymbol{u} - \varepsilon^p)] \, dV \\ & - \int_{\Omega_p} \delta \boldsymbol{\sigma} : [\varepsilon^p - \varepsilon_n^p - \lambda \partial_\sigma \phi] \, dV - \int_{\Omega_p} \delta \lambda \, dV \\ & + \int_{\Omega_p} \delta \boldsymbol{q} \cdot [-\boldsymbol{\alpha} + \boldsymbol{\alpha}_n - \lambda \partial_\boldsymbol{q} \phi] \, dV = 0 \end{aligned} \quad (6) $$

holds true for any displacement variation $\delta \boldsymbol{u}$ solving the homogeneous elastic equilibrium field equation, i.e. of the form

$$ \delta u_k(\boldsymbol{y}) = \int_{S_u} \delta \tilde{t}_i(\boldsymbol{x}) U_i^k(\boldsymbol{x}-\boldsymbol{y}) \, dS_x - \int_{S_t} \sum_{ij}^k (\boldsymbol{x}-\boldsymbol{y}) n_j(\boldsymbol{x}) \delta \tilde{u}_i(\boldsymbol{x}) \, dS_x \quad (7) $$
(with $\delta\tilde{u} \in \tilde{H}^{1/2}(S_t)$ and $\delta\tilde{t} \in H^{-1/2}(S_u)$), stress variations $\delta\sigma \in L^2(\Omega)$, plastic strain variations $\delta\varepsilon^p \in L^2(\Omega)$ and variations $\delta\alpha \in L^2(\Omega)$, $\delta q \in L^2(\Omega)$, $\delta\lambda \in \{\gamma \in L^2(\Omega), \gamma \ge 0\}$. From eq. (6) on, $(\tilde{u}, \tilde{t}, \sigma, \varepsilon^p, q, \alpha)$ are implicitly taken at $t = t_{n+1}$ and the index $n+1$ is omitted for brevity.

The substitution of (7) in (6) mandates a careful treatment, because the traces of $\delta u$ and $[C : \nabla^s \delta u] \cdot n$ on $S$ involve non-integrable kernels. To overcome this difficulty, the following definition of the hypersingular boundary operator on $S$ is used [8]

$$ \int_S \phi \cdot (C : \nabla^s \delta u) \cdot n \, dS = \int_{\Omega} \nabla^s (\gamma_0^- \phi) : C : \nabla^s \delta u \, dV \quad \forall \phi \in L^2(S) \quad (8) $$

for $\delta u$ satisfying $\mathrm{div}(C : \nabla^s \delta u) = 0$ in $\Omega$ ($\gamma_0^- : \mathrm{extension \ operator \ around} S$). By taking into account the following decomposition of the hypersingular kernel $C_{ijab}\Sigma_{kl,b}^a(r)$ [13]

$$ C_{ijab}\Sigma_{kl,b}^a(r) = -e_{jfq}e_{lhs}B_{iqks,fh}(r) - C_{ijkl}\delta(0) \\ B_{iqks}(r) = \frac{G}{4\pi(1-\nu)r} [\delta_{qs}r,i r,k + \delta_{ik}\delta_{qs} - 2\nu\delta_{is}\delta_{kq} - (1-\nu)\delta_{iq}\delta_{ks}] $$

and applying Stokes theorem to reduce the singularity of kernels in both boundary and domain integral operators, a fully regularized variational formulation is finally obtained from (6):

$$ \begin{align} a_{uu}(\delta \tilde{u}, \tilde{u}) + a_{ut}(\delta \tilde{u}, \tilde{t}) &= \ell_u(\delta \tilde{u}) + a_{u\epsilon}(\delta \tilde{u}, \epsilon^p) \nonumber \\ a_{tu}(\delta \tilde{t}, \tilde{u}) + a_{tt}(\delta \tilde{t}, \tilde{t}) &= \ell_t(\delta \tilde{t}) + a_{t\epsilon}(\delta \tilde{t}, \epsilon^p) \nonumber \\ \int_{\Omega_p} \delta \epsilon^p : \sigma \, dV + a_{\epsilon\epsilon}^{\sigma}(\delta \epsilon^p, \epsilon^p) &= \ell_{\epsilon}(\delta \epsilon^p) + a_{\epsilon u}(\delta \epsilon^p, \tilde{u}) + a_{\epsilon t}(\delta \epsilon^p, \tilde{t}) \tag{9} \\ &\quad + \text{Evolution laws in weak form} \nonumber \end{align} $$

with some properties of symmetry and sign-definiteness proved in [16]. The terms corresponding to the purely elastic response of the body can
be found in [2] or [10], while the new ones are given by:

$$a_{u\epsilon}(\delta\tilde{u}, \epsilon^p) = - \int_{S_t} \int_{\Omega_p} R_s \delta\tilde{u}_k(\mathbf{x}) e_{jfq} B_{iqks,f}(r)(\mathbf{y}-\mathbf{x}) \epsilon_{ij}^p(\mathbf{y}) \, dS_y \, dS_x$$

$$a_{t\epsilon}(\delta\tilde{t}, \epsilon^p) = \int_{S_u} \int_{\Omega_p} \delta\tilde{t}_k(\mathbf{x}) \Sigma_{ij}^k(\mathbf{y}-\mathbf{x}) \epsilon_{ij}^p(\mathbf{y}) \, dV_y \, dS_x$$

$$\begin{aligned} a_{\epsilon\epsilon}^{\sigma}(\delta\epsilon^p, \epsilon^p) = & \int_{\Omega_p} \int_{\Omega_p} (e_{qjf} \delta\epsilon_{ij,f}^p)(\mathbf{x}) B_{iqks}(\mathbf{y}-\mathbf{x}) (e_{shl} \epsilon_{kl,h}^p)(\mathbf{y}) \, dV_y \, dV_x \\ & - \int_{\Omega_p} \int_{S_p} (e_{qjf} \delta\epsilon_{ij,f}^p)(\mathbf{x}) B_{iqks}(\mathbf{y}-\mathbf{x}) (e_{shl} n_h \epsilon_{kl}^p)(\mathbf{y}) \, dS_y \, dV_x \\ & - \int_{\Omega_p} \int_{S_p} (e_{qjf} \epsilon_{ij,f}^p)(\mathbf{x}) B_{iqks}(\mathbf{y}-\mathbf{x}) (e_{shl} n_h \delta\epsilon_{kl}^p)(\mathbf{y}) \, dS_y \, dV_x \\ & + \int_{S_p} \int_{S_p} (e_{qjf} \delta\epsilon_{ij}^p n_f)(\mathbf{x}) B_{iqks}(\mathbf{y}-\mathbf{x}) (e_{shl} \epsilon_{kl}^p n_h)(\mathbf{y}) \, dS_y \, dS_x \end{aligned}$$

Following [16], this set of variational equations is found to express the stationarity of the functional $L(\chi^\sigma)$ defined by

$$\begin{aligned} L(\chi^\sigma) = & \frac{1}{2} a_{uu}(\tilde{u}, \tilde{u}) + \frac{1}{2} a_{tt}(\tilde{t}, \tilde{t}) + \frac{1}{2} a_{\epsilon\epsilon}^{\sigma}(\epsilon^p, \epsilon^p) \\ & + a_{ut}(\tilde{u}, \tilde{t}) - a_{u\epsilon}(\tilde{u}, \epsilon^p) - a_{t\epsilon}(\tilde{t}, \epsilon^p) - \ell_u(\tilde{u}) - \ell_t(\tilde{t}) - \ell_{\epsilon}(\epsilon^p) \\ & + \int_{\Omega_p} \sigma : (\epsilon^p - \epsilon_n^p) \, dV + \int_{\Omega_p} [ \Theta(\alpha) - q.(\alpha - \alpha_n) ] \, dV \\ & - \int_{\Omega_p} \lambda \, \phi(\sigma, q) \, dV \end{aligned}$$

Note that $L(\chi^\sigma)$ does not involve the total strain $\epsilon$, so that the present formulation might be termed an *assumed-stress method*. Some important results about consistency, stability and convergence of the associated numerical scheme has been established in [7].

Performing the change of variable $\sigma = C : (\epsilon - \epsilon^p)$ would yield an *assumed-strain method*, commonly used in FEM and also introduced for collocation BEM in [4]. Here, this manipulation changes $L(\chi^\sigma)$ into a new functional $L(\chi^\epsilon)$ where $\chi^\epsilon = ((\tilde{u}, \tilde{t}), \epsilon, \epsilon^p, q, \alpha)$. The resulting assumed-strain formulation, which will be retained in the following sections, has obviously the same properties as the assumed-stress one.

### 3. LOCAL INTEGRATION

The inelastic behaviour of the body is here assumed to be described by a Von-Mises yield criterion associated to a linear isotropic and kinematic
hardening rule respectively characterized by two positive constants $h$ and $H$ (see [21] for more general classes of materials). The evolution laws resulting of the stationarity conditions of $L(\chi^\epsilon)$ are written in a weak form as:

$$ 
\begin{aligned}
\int_{\Omega_p} \delta\epsilon : \mathbf{C} : (\varepsilon^p - \varepsilon_n^p) \, dV &= \int_{\Omega_p} \delta\epsilon : \mathbf{C} : \lambda \hat{n} \, dV \\
\int_{\Omega_p} \delta\beta : (\beta - \beta_n) \, dV &= \int_{\Omega_p} \delta\beta : \frac{2H}{3} \lambda \hat{n} \, dV \\
\int_{\Omega_p} \delta p : (p - p_n) \, dV &= \int_{\Omega_p} \delta p : \sqrt{\frac{2}{3}} \lambda \, dV \\
\int_{\Omega_p} \delta\lambda\phi \, dV &= 0
\end{aligned}
\quad (10) $$

with the classical notations

$$ 
\begin{aligned}
\phi(\sigma - \beta, p) &= \| s - \beta \| - \sqrt{\frac{2}{3}(\sigma_0 + hp)} \\
\hat{n} &= \frac{\xi}{\| \xi \|} \quad \text{with} \quad \xi = s - \beta \\
s &= s_n + 2G(e - e_n - \lambda \hat{n})
\end{aligned}
\quad (11) $$

where $\beta$ is the center of the elastic domain, $p$ the cumulated plastic strain, $s$ and $e$ the deviatoric stress strain, respectively.

These set of relations are solved by considering the strain-driven Return Mapping Algorithm (RMA), widely used in FEM [21]. The RMA is based on a predictor-corrector scheme and allows to construct an intermediate configuration consistent with the normality rules and the Hooke's law. If $\lambda = 0$ the material remains elastic, otherwise the internal variables must be updated. $\lambda$ results from the consistency condition at $t_{n+1}$, which corresponds to the projection onto the yield surface of a trial stress defined by $\xi^E$, and is given by:

$$ 
\begin{aligned}
\xi^E &= s_n + 2G(e - e_n) - \beta_n, \quad \hat{n} = \frac{\xi^E}{\| \xi^E \|} \quad \text{and} \\
\lambda &= \frac{3}{2(3G + H + h)} \langle \phi(\xi^E, p_n) \rangle
\end{aligned}
\quad (12) $$

A similar scheme corresponding to a stress-driven RMA is described in [19] but the plastic multiplier $\lambda$ can not be evaluated in the case of perfect plasticity.
4. NUMERICAL IMPLEMENTATION

The numerical solution procedure is based upon a discretization of eqs. (9) and (10) in space. As pointed out in [15], the usual boundary element interpolation of unknown displacements and tractions on the boundary (which model the elastic structural behaviour) must be supplemented with a domain interpolation of strains (the potentially plastic region being cut into ‘cells’) when non-linear material behaviour is present. Similar considerations arise in collocation BEM as well [4, 6].

$$
\begin{align}
& y(\boldsymbol{x}) = \boldsymbol{N}(\boldsymbol{x}).\boldsymbol{Y}, \quad \boldsymbol{\varepsilon}(\boldsymbol{x}) = \boldsymbol{N}(\boldsymbol{x}).\boldsymbol{E}, \quad \boldsymbol{\varepsilon}^p(\boldsymbol{x}) = \boldsymbol{N}(\boldsymbol{x}).\boldsymbol{E}^p \nonumber \\
& b(\boldsymbol{x}) = \boldsymbol{N}(\boldsymbol{x}).\boldsymbol{B}, \quad p(\boldsymbol{x}) = \boldsymbol{N}(\boldsymbol{x}).\boldsymbol{P} \tag{13} \\
& \text{and } \lambda(\boldsymbol{x}) = \max(\boldsymbol{N}(\boldsymbol{x}).\boldsymbol{\Lambda}, 0) \nonumber
\end{align}
$$

where (E, E<sup>p</sup>, B, P, Λ) denotes the vector of domain nodal unknowns and Y the vector of boundary nodal densities. In this way we consider a continuous interpolation of the unknown boundary densities y = (tilde{u}, tilde{t}), the total strain ε, the plastic strain ε and the internal variables β and p, contrary to [16] or [11]. In the last cited paper, the authors used compatible plastic strain field in order to avoid domain integral calculation.

A Galerkin scheme is then applied to equations (9) and to the time-
discretized plasic strain flow rule. The evolutions of β, p and the consis-
tency condition are performed at the *nodes* of the domain mesh, which
amounts to use Dirac measures at the nodes of Ω<sub>p</sub> for the test functions
in eqs. (10). The resulting non-linear algebraic system is written as

$$
K.Y_{n+1} = F_{n+1} + Q.E_{n+1}^{p} \quad (14)
$$

$$
M.E_{n+1} = H_{n+1} + Q^t.Y_{n+1} - Z.E_{n+1}^p \quad (15)
$$

$$
M.E_{n+1}^{p} = M.E_{n}^{p} + \Psi_{\varepsilon}(E_{n+1}^{p}, E_{n}, E_{n}^{p}) \quad (16)
$$

+ Evolution of internal variables and consistency condition

with the notations

$$
\Psi_{\varepsilon}(\mathbf{E}_{n+1}) = 2G \int_{\Omega_p} \mathbf{N}^t(\mathbf{x}) \lambda_{n+1}(\mathbf{x}) \hat{\mathbf{n}}_{n+1} \, \mathrm{d}V \quad (17)
$$

$$
M = \int_{\Omega_p} \mathbf{N}^t(\mathbf{x}) : \mathbf{C} : \mathbf{N}(\mathbf{x}) \, \mathrm{d}V \quad (18)
$$

Solving eq. (14) for $\mathbf{Y}_{n+1}$ and substituting the result into eqs. (15,16) yields a non-linear system of equations for $\mathbf{E}_{n+1}$ and $\mathbf{E}_{n+1}^p$ characterizing the material plastic constitutive behaviour:

$$
M.E_{n+1} = G_{n+1} - \Xi.E_{n+1}^{p} \quad (19)
$$

$$
M.E_{n+1}^{p} = M.E_{n}^{p} + \Psi_{\varepsilon}(E_{n+1}, E_{n}, E_{n}^{p}) \quad (20)
$$
which formally involves a nonlinear vector function $\mathcal{A}_\epsilon$ whose tangent operator is the global consistent tangent operator (CTO) [20]:

$$ \mathcal{A}_{\epsilon} \begin{pmatrix} \mathbf{E}_{n+1}^{p} \\ \mathbf{E}_{n+1} \end{pmatrix} = \begin{pmatrix} \mathbf{u}_{n+1}^{0} \\ \mathbf{t}_{n+1}^{0} \end{pmatrix} \qquad \partial \mathcal{A}_{\epsilon} = \begin{bmatrix} \Xi & \mathbf{M} \\ \mathbf{M} & \Gamma(\mathbf{E}_{n+1}) \end{bmatrix} \\ \Gamma(\mathbf{E}_{n+1}) \equiv \mathbf{M} \cdot \frac{\partial \mathbf{E}_{n+1}^{p}}{\partial \mathbf{E}_{n+1}} = 2G \int_{\Omega_p} \mathbf{N}^{t}(\mathbf{x}) \frac{\partial(\lambda_{n+1}\hat{\mathbf{n}}_{n+1})}{\partial \varepsilon_{n+1}} \mathbf{N}(\mathbf{x}) \, dV $$

The global CTO $\partial \mathcal{A}_\epsilon$ is symmetric but from a numerical point of view, this formulation leads to an algebraic system of large size. Unfortunately neither $\Gamma$ nor $\Xi$ are invertible [5, 14]. It is more interesting to eliminate the plastic strain vector $\mathbf{E}_{n+1}^p$ in order to get a non-linear equation for $\mathbf{E}_{n+1}$. Equations (19,20) then become (with $\mathbf{E} \equiv \mathbf{E}_{n+1}$):

$$ \mathbf{R}_{\epsilon}(\mathbf{E}) \equiv \mathbf{E} + \mathbf{M}^{-1}\Xi\mathbf{M}^{-1}(\Psi_{\epsilon}(\mathbf{E}) + \mathbf{M}\mathbf{E}_{n}^{p}) - \mathbf{M}^{-1}\mathbf{G}_{n+1} = 0 \quad (21) $$

This equation is solved with a Newton-Raphson method and the associated consistent iterative scheme reads:

$$ [\mathbf{I} + 2\mathbf{G}\mathbf{M}^{-1}\Xi\mathbf{M}^{-1}\Gamma]\Delta\mathbf{E}^i = -\mathbf{R}_\epsilon(\mathbf{E}^i) \qquad \mathbf{E}^{i+1} = \mathbf{E}^i + \Delta\mathbf{E}^i \quad (22) $$

where $\mathbf{I} + 2\mathbf{G}\mathbf{M}^{-1}\Xi\mathbf{M}^{-1}\Gamma$, the global CTO corresponding to (22), is seen to be non-symmetric and fully populated because of $\Xi$. This presentation encompasses the treatment proposed in [11] which corresponds to a modified Newton-Raphson scheme with a tangent operator equal to $\mathbf{I}$.

Concerning the numerical evaluation of the integral operators, since the present formulation involves integrable kernels, classical Gauss quadrature rules are used for the inner integration (after removing of the singularity by Duffy's coordinates) and the outer integration (the inner potential being regular). This method is easy to implement but numerical quadrature errors entail slight loss of symmetry because the inner and outer integrations are not treated symmetrically [17].

## 5. NUMERICAL RESULTS

As a test example, the elastoplastic torsion problem for a cylindrical body (radius $r = 100$ mm, height $h = 50$ mm) is considered. A displacement-controlled twisting motion is prescribed on the plane ends. Perfect plasticity is considered; the material properties are $E = 200000$ MPa, $\nu = 0.3$ and $\sigma_0 = 240$ MPa. The numerical model is made of 8 27-noded hexagonal cells and 48 9-noded quadrilateral boundary elements (Fig. 1).

The results are compared to the analytical solution and the load-displacement curve is plotted in Fig. 2.
<figure><img src="image_2.png" /><figcaption>Figure 1 Discretized BE-CE model. Boundary and internal mesh</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2 Twisting moment-Angular rotation curve</figcaption></figure>

## 6. CONCLUSION

A BEM-based strain formulation to quasi-static elastoplasticity evolutive analysis has been presented for general standard materials. The associated iterative scheme admits a global symmetric CTO. Some numerical aspects have now to be improved: quadrature rules taking into account the symmetry of the integration domains and parallel computing. Extensions to damage or finite strain models can be performed in the same way.
References

[1] M. Bonnet. *Boundary integral equation methods in solids and fluids*. Wiley, 1999.

[2] M. Bonnet. Regularized direct and indirect symmetric variational bie formula-
tions for three-dimensional elasticity. *Eng. Anal. with Bound. Elem.*, 15:93-102,
1995.

[3] M. Bonnet, G. Maier, and C. Polizzotto. Symmetric galerkin boundary element method. *Appl. Mech. Review*, 51:669-704, 1998.

[4] M. Bonnet and S. Mukherjee. Implicit bem formulations for usual and sensitivity problems in elastoplasticity using the consistent tangent operator concept. *Int. J. Solids Struct.*, 33:4461-4480, 1996.

[5] H.D. Bui and K. Dangvan. Sur le problème aux limites en vitesse des contraintes du solide élasto-plastique. *Int. J. Solids Struct.*, 6:183-193, 1970.

[6] B. Burgardt and A. Levan. A fully regularized boundary formulation for three-dimensional elastoplastic problems. In C. Brebbia, A. Kassab, and M. Choppa, editors, *Boundary Elements 20*, pages 95-103. Computational Mechanics Publications, 1998.

[7] C. Comi and G. Maier. Extremum, convergence and stability properties of the finite-increment problem in elastic-plastic boundary element method. *Int. J. Solids Struct.*, 29:249-270, 1992.

[8] M. Costabel. Boundary integral operators on lipschitz domains : elementary results. *SIAM J. Math. Anal.*, 19:613-626, 1988.

[9] J. Lemaitre and J.C. Chaboche. *Mechanics of solid materials*. Cambridge, 1994.

[10] S. Li, M. E. Mear, and L. Xiao. Symmetric weak-form integral equation method for three-dimensional fracture analysis. *Comput. Meth. in Appl. Mech. Eng.*, 51:669-704, 1998.

[11] G. Maier, S. Miccoli, G. Novati, and U. Perego. Symmetric galerkin boundary element method in plasticity and gradient plasticity. *Comput. Mech.*, 17:115-129, 1995.

[12] G. Maier and C. Polizzotto. A galerkin approach to boundary element elasto-plastic analysis. *Comput. Meth. in Appl. Mech. Eng.*, 60:175-194, 1987.

[13] J.C. Nedelec. Integral equations with non-integrable kernels. *Integral equations and operator theory*, 5:562-572, 1982.

[14] C. Polizzotto. An energy approach to the boundary element method. part ii : Elastic-plastics solids. *Comput. Meth. in Appl. Mech. Eng.*, 69:263-276, 1988.

[15] C. Polizzotto, C. Panzeca, and M. Zito. A consistent boundary/interior element method for evolutive elastic-plastic structural analysis. In J.H. Kane, G. Maier, N. Tosaka, and S.N. Atluri, editors, *Advances in Boundary Element Techniques*, pages 349-369. Springer Verlag, 1993.

[16] C. Polizzotto and M. Zito. A step-wise variational approach to elastic-plastic analysis by boundary-interior elements. *Int. J. of Plasticity*, 10:81-101, 1994.

[17] S. A. Sauter and C. Lage. Transformation of hypersingular integrals and black-box cubature. Technical Report 97-17, Universität Kiel, 1997. Extended version.
[18] J. C. Simo and S. Govindjee. Non-linear b-stability and symmetry preserving return mapping algorithms for plasticity and viscoplasticity. *Int. J. for Num. Meth. in Engng.*, 31:151–176, 1991.

[19] J. C. Simo, J. G. Kennedy, and R. L. Taylor. Complementary mixed finite element formulations for elastoplasticity. *Comput. Meth. in Appl. Mech. Engng.*, 74:177–206, 1989.

[20] J. C. Simo and R. L. Taylor. Consistent tangent operators for rate independent elastoplasticity. *Comput. Meth. in Appl. Mech. Engng.*, 48:101–118, 1985.

[21] J.C. Simo and T.J.R. Hughes. *Computational Inelasticity*, volume 7 of *Interdisciplinary Applied Mathematics*. Springer Verlag, 1998.