# Accounting for Incompressibility in Reproducing Kernel Particles Meshless Approximations

Pierre Joyot, Jean Trunzler, Francisco Chinesta

► To cite this version:

Pierre Joyot, Jean Trunzler, Francisco Chinesta. Accounting for Incompressibility in Reproducing Kernel Particles Meshless Approximations. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2006, 15 (5), pp.513-527. 10.3166/remn.15.513-527. hal-01004971

HAL Id: hal-01004971

https://hal.science/hal-01004971

Submitted on 26 Mar 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.

Public Domain
Accounting for incompressibility
in reproducing kernel particle meshless
approximations

Pierre Joyot* — Jean Trunzler* — Francisco Chinesta**

* LIPSI-ESTIA

Technopole Izarbel, F-64210 Bidart

{j.trunzler, p.joyot}@estia.fr

** LMSP UMR 8106 CNRS-ENSAM

151 Boulevard de l'Hôpital, F-75013 Paris

francisco.chinesta@paris.ensam.fr

**ABSTRACT.** Meshless approximations seem to be an appealing choice for simulating forming processes involving large transformations because they allows alleviating the mesh constraints. However, because the novelty of these techniques a lot of questions are today unresolved. One of these open problems is the treatment of incompressibility which as well known impose some restrictions on the choice of the approximation spaces. The accurate treatment of incompressibility is a key point in the simulation of forming processes because the plastic flow can be in fact considered as incompressible. This paper introduces the problematic as well as some possibilities for taking into account the incompressibility in the context of mixed formulations, making special emphasis in a kind of Hermite approximations.

**RÉSUMÉ.** Les méthodes sans maillage semblent bien adaptées à la résolution de problèmes en grandes transformations, mais à cause de leur jeunesse quelques difficultés subsistent. Une de ces difficultés concerne la prise en compte de l'incompressibilité qui induit des contraintes sur le choix des espaces fonctionnels d'approximation. La prise en compte de l'incompressibilité est fondamentale dans le cadre de la simulation des procédés de mise en forme. Dans ce travail nous présentons la problématique liée à l'incompressibilité matérielle ainsi que différentes possibilités pour sa prise en compte dans le cadre des formulations mixtes, où des approximations de type Hermite sont particulièrement adaptées.

**KEYWORDS:** RKPA, meshless approximation, mixed formulations, incompressibility.

**MOTS-CLÉS :** RKPA, méthodes sans maillage, approximations mixtes, incompressibilité.
# 1. Introduction

Meshless methods are an appealing choice for constructing functional approximations (with different degrees of consistency and continuity) without a mesh support. Thus, this kind of techniques seem to be specially appropriated for treating 3D problems involving large displacements, due to the fact that the approximation is constructed only from the cloud of nodes whose positions evolve during the material deformation. In this manner neither remeshing nor fields projections are a priori required.

Other important point lies in the easy introduction of some known information related to the problem solution within the approximation functional basis. For this purpose, different reproduction conditions are enforced in the construction of the approximation functions. This approach has been widely used in the context of the moving least squares approximations involved in the diffuse meshless techniques (Nayroles et al., 1992) as well as in the element free Galerkin method (Belytschko et al., 1996).

Meshless approximations seem to be an appealing choice for simulating forming processes involving large transformations because they allows alleviating the mesh constraints. However, because the novelty of these techniques a lot of questions are today unresolved. One of these open problems is the treatment of incompressibility which as well known impose some restrictions on the choice of the approximation spaces. The accurate treatment of incompressibility is a key point in the simulation of forming processes because the plastic flow can be in fact considered as incompressible. This paper introduces the problematic as well as some possibilities for taking into account the incompressibility in the context of mixed formulations, making special emphasis in a kind of Hermite approximations.

# 2. Reproducing Kernel Particle Approximation

The Reproducing Kernel Particle Method is based in the following approximation:

$$f^a(x) = \int_{\Omega} \phi(x-s, h) f(s) ds \quad [1]$$

where h defines the size of the support of the kernel function, which approaches the delta distribution as $h \to 0$.

If one want to reproduce polynomials up to a certain degree the kernel function must be constructed in an appropriate manner. For the sake of simplicity we consider the 1D case and second degree reproducing conditions, which imply:
$$ egin{cases} 1 = \int_{\Omega} \phi(x-s,h)1ds \\ x = \int_{\Omega} \phi(x-s,h)sds \\ x^2 = \int_{\Omega} \phi(x-s,h)s^2ds \end{cases} \quad [2] $$

As usual kernel functions do not verify the reproducing conditions [2] (Liu et al., 1995) proposed to introduce a correction in the form:

$$ f^a(x) = \int_{\Omega} C(x, x-s)\phi(x-s, h)f(s)ds \quad [3] $$

with the correction function given by:

$$ C(x, x-s) = \underline{H}^T(x-s)\underline{b}(x) \quad [4] $$

where

$$ \underline{H}^T(x-s) = [1, (x-s), (x-s)^2] \quad [5] $$

and the components of $\underline{b}(x)$ must be determined for verifying the reproduction conditions.

The discrete form of Equation [3] results:

$$ f^a(x) = \sum_J C(x, x-x_J)\phi(x-x_J, h)f(x_J) \quad [6] $$

with

$$ C(x, x-x_J) = \underline{H}^T(x-x_J)\underline{b}(x) \quad [7] $$

Thus, the reproduction conditions result:
$$ egin{cases} 1 = \sum_J \underline{H}^T (x - x_J) \underline{b}(x) \phi(x - x_J, h) \\ x = \sum_J \underline{H}^T (x - x_J) \underline{b}(x) \phi(x - x_J, h) x_J \\ x^2 = \sum_J \underline{H}^T (x - x_J) \underline{b}(x) \phi(x - x_J, h) x_J^2 \end{cases} \quad [8] $$

that can be written in the following matrix form:

$$ \underline{R}(x) = \left[ \sum_J \underline{R}(x_J) \underline{H}^T (x - x_J) \phi(x - x_J, h) \right] \underline{b}(x) \quad [9] $$

or

$$ \underline{R}(x) = \underline{M}(x) \underline{b}(x) \quad [10] $$

where $\underline{R}$ denotes the reproduction vector and $\underline{M}(x)$ the moment matrix. Thus, it results:

$$ \underline{b}(x) = \underline{M}^{-1}(x) \underline{R}(x) \quad [11] $$

and

$$ \begin{aligned} f^a(x) &= \sum_J \underbrace{\underline{H}^T (x - x_J) \underline{M}^{-1}(x) \underline{R}(x)}_{C(x, x-x_J)} \phi(x - x_J, h) f(x_J) = \\ &= \sum_J N_J(x) f(x_J) \end{aligned} \quad [12] $$

where $N_J(x)$ denote the approximation shape functions that do not verify the delta property.

**REMARK 1.** — This procedure can be generalized for reproducing arbitrary functions related to the problem solution (Trunzler *et al.*, 2005).

**REMARK 2.** — In the 2D case the reproduction vector results

$$ \underline{R} = [1, x, y, x^2, xy, y^2] \quad \text{for the same degree of reproduction.} $$

**REMARK 3.** — The integration weights are not present in Equation [6] because they are implicitly considered in the expression of $\underline{b}(x)$.
REMARK 4. — The trial and test functions related to the variational formulation of a generic problem can be approximated using the scheme [6] from which the discrete problem is derived as in the usual finite element method. However due to the non polynomial character of the kernel function richer integration rules must be considered.

The main drawbacks related to this kind of discretization technique is the difficulty to impose essential boundary condition (the shape functions related to internal nodes don't vanish on the domain boundary) as well as the difficulty to control the shape function support when the nodal distribution evolves significantly (as encountered for example in forming processes simulations).

Moreover, when this kind of approximation is considered for treating incompressible media one must verify the LBB stability conditions which implies, when a mixed variational formulation is considered, richer approximation in the velocity or displacement fields than in the pressure one (which is in fact the Lagrange multiplier related to the incompressibility restriction).

In the context of our former works, we have proved that by changing the reproduction condition (second degree for the richer field and first degree for the one related to the Lagrange multiplier) the derived approximation is not stable. The same conclusion was obtained by changing simultaneously the shape functions support. Specific strategies have been proposed for solving some kind of problems, as the Stokes one for example, but its generality is limited. For this reason we introduce in this paper richer approximations based on the use of Hermite-RKPA approximations.

### 3. Efficient Hermite – RKPM approximations

In order to derive a more general H-RKPA approximation, we consider (in the 1D case for the sake of simplicity):

$$ \left( \frac{f^a(x)}{df^a(x)} \right) = \sum_{J=1}^{N} C \phi \left( \frac{f(x_J)}{df(x_J)} \right) \quad [13] $$

If we are looking for second degree reproduction condition, we can write:

$$ \left( \begin{array}{c} a + bx + cx^2 \\ b + 2cx \end{array} \right) = \sum_{J=1}^{N} C \phi \left( \begin{array}{c} a + bx_J + cx_J^2 \\ b + 2cx_J \end{array} \right) \quad [14] $$

that implies:
$$ \begin{pmatrix} 1 & x & x^2 \\ 0 & 1 & 2x \end{pmatrix} = \sum_{J}^{N} \underline{\underline{C}} \phi \begin{pmatrix} 1 & x_J & x_J^2 \\ 0 & 1 & 2x_J \end{pmatrix} \quad [15] $$

or

$$ \begin{pmatrix} 1 & 0 \\ x & 1 \\ x^2 & 2x \end{pmatrix} = \sum_{J}^{N} \underline{\underline{\phi}} \begin{pmatrix} 1 & 0 \\ x_J & 1 \\ x_J^2 & 2x_J \end{pmatrix} \underline{\underline{C}}^T \quad [16] $$

Defining now the reproduction matrix

$$ \underline{\underline{R}}(x) = \begin{pmatrix} 1 & 0 \\ x & 1 \\ x^2 & 2x \end{pmatrix} \quad [17] $$

the correction matrix is expressed by:

$$ \underline{\underline{C}}^T = \underline{\underline{H}}^T \underline{\underline{b}} \quad [18] $$

with

$$ \underline{\underline{H}}(x - x_J) = \begin{pmatrix} 1 & 0 \\ x - x_J & 1 \\ (x - x_J)^2 & 2(x - x_J) \end{pmatrix} \quad [19] $$

Thus, from Equations [16] and [18] it results:

$$ \underline{\underline{R}}(x) = \left[ \sum_{J}^{N} \underline{\underline{\phi}} \underline{\underline{R}}(x_J) \underline{\underline{H}}^T \right] \underline{\underline{b}} = \underline{\underline{M}} \underline{\underline{b}} \quad [20] $$

which allows to compute $\underline{\underline{b}}, \underline{\underline{b}} = \underline{\underline{M}}^{-1} \underline{\underline{R}}(x)$. Thus, we can finally write:

$$ \left( \frac{f^a(x)}{df^a(x)} \right) = \sum_{J}^{N} \underline{\underline{b}}^T \underline{\underline{H}} \underline{\underline{\phi}} \left( \frac{f(x_J)}{df(x_J)} \right) = \sum_{J}^{N} \underline{\underline{N}} \left( \frac{f(x_J)}{df(x_J)} \right) \quad [21] $$
REMARK 5. — It is easy to verify that the derivatives involved in the previous expression are diffuse derivatives that can be employed in collocation discretization schemes but not discretization schemes operating on the variational formulation. However, usual collocation schemes using the nodes as collocation points lead to a rectangular discrete linear system, because there are several dof at each node.

# 4. Numerical results

## 4.1. Approximation tests

In this section we analyze the approximation of the polynomials. Firstly, we consider an approximation with third degree of consistency, i.e.

$$ \underline{R}(x) = \begin{pmatrix} 1 & 0 \\ x & 1 \\ x^2 & 2x \\ x^3 & 3x^2 \end{pmatrix} $$

involving two degrees of freedom at each node, one related to the unknown field and the other one to its derivative. When the approximated polynomial has a degree lower or equal to the one considered in the reproduction conditions, that polynomial is approximated with the machine precision. Figure 1 depicts the approximated (stars) versus the exact (solid curve) solutions related to the polynomial:

$$ f^{\text{exact}}(x) = x - 5x^2 + 5x^3 $$

<figure><img src="image_1.png" /><figcaption>Figure 1. Approximation of polynomial $f(x) = x - 5x^2 + 5x^3$</figcaption></figure>
Figure 2 depicts similar results concerning the derivative.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Approximation of the derivative of $f(x) = x - 5x^2 + 5x^3$</figcaption></figure>

As the polynomial and its derivative are exactly reproduced, the diffuse derivative corresponds to the exact one, which justifies the excellent accuracy obtained in the derivatives. Errors are of the order of the machine precision.

Finally, we are analyzing the convergence of the approximation when a polynomial that cannot be exactly reproduced: $f(x) = -x^4 + x^5$. Figure 3 depicts the computed solution versus the exact one.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Approximation of polynomial $f(x) = -x^4 + x^5$</figcaption></figure>
If we compute the approximation error using the L² norm, we obtain fourth order convergence with respect to the log of the total number of degrees of freedom.

## 4.2. Second order PDE

In order to analyze the behavior of the H-RKPM approximation just defined when it is used for solving second order partial differential equation, we first consider the 1D problem defined in [0,1]:

$$ \frac{d^2 f}{dx^2} = -12x^2 + 6x \text{ with } f(x=0) = 0 \text{ and } f(x=1) = 0 $$

whose exact solution results $f(x) = -x^4 + x^5$.

The approximation is built from the reproduction matrix:

$$ \underline{R}(x) = \begin{pmatrix} 1 & 0 \\ x & 1 \\ x^2 & 2x \end{pmatrix} $$

Figure 4 depicts the exact and the computed solution when 20 nodes are considered, each one involving 2 degrees of freedom: the field and its derivative.

<figure><img src="image_4.png" /><figcaption>Figure 4. Computed (stars) versus exact (solid curve) solution</figcaption></figure>
The convergence order obtained was only first order, significantly lower than the third order expected when the L² norm is considered, may be induced by the numerical integration.

REMARK 6. — The boundary conditions have been imposed by removing the rows of the linear system related to the first degree of freedom of the first and last nodes, imposing at its place the Dirichlet boundary conditions. If these conditions are imposed by using Lagrange multipliers the solution improvement is inappreciable.

REMARK 7. — The variational formulation can be integrated with or without integration by parts. Obviously, due to the approximation continuity and derivability no integration by parts is required. In the case of integration by parts the boundary integrals must be considered because the test functions do not vanish on the domain boundary. In both cases the solution accuracy does not exhibit significant differences.

REMARK 8. — When a regular distribution of quadrature points is used for integrating the variational formulation the shape function derivatives can be computed from the value of the shape functions at the integration points using finite differences. In practical applications the accuracy is not significantly degraded.

REMARK 9. — When the variational formulation is considered, as argued previously, we can not use for defining the field derivative its diffuse derivative given by the second row of Equation [21].

REMARK 10. — This technique can be applied directly for solving fourth order PDE which requires the extra imposition of some field derivative on the domain boundary with either collocation or Lagrange multipliers.

We consider now the 2D problem defined in [0,1]²:

$$ \Delta f = 4 \quad \text{with} \quad f(x=0, y) = y^2, \quad f(x, y=0) = x^2, \quad f(x=1, y) = 1+y^2 \quad \text{and} \quad f(x, y=1) = 1+x^2, \quad \text{whose exact solution is given by} \quad f(x, y) = x^2 + y^2. $$

In this case the reproduction matrix results

$$ \underline{R}(x, y) = \begin{pmatrix} 1 & 0 & 0 \\ x & 1 & 0 \\ y & 0 & 1 \\ x^2 & 2x & 0 \\ xy & y & x \\ y^2 & 0 & 2y \end{pmatrix} $$

which has three nodal degrees of freedom: the unknown field and its first derivatives with respect to both space coordinates. The Dirichlet conditions are imposed, as previously, *via* collocation. Figure 5 depicts the numerical versus de computed
solutions. As both solutions coincide (at the machine precision) they are indistinguishable in the representation.

<figure><img src="image_5.png" /><figcaption>Figure 5. Computed and exact solutions: both are perfectly superposed</figcaption></figure>

### 4.3. Accounting for incompressibility: the Stokes problem

This approximation can be applied for stabilizing the discretization of the Stokes problem in $[0,1]^2$:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial p}{\partial x} + \Delta u = 0 \\ \\ \displaystyle \frac{\partial p}{\partial y} + \Delta v = 0 \\ \\ \displaystyle \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \end{array} \right. \qquad [22] $$

whose variational formulation (without integration by parts) results:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \left( \frac{\partial p}{\partial x} + \Delta u \right) u^* = 0 \\ \\ \displaystyle \int_{\Omega} \left( \frac{\partial p}{\partial y} + \Delta v \right) v^* = 0 \\ \\ \displaystyle \int_{\Omega} \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) p^* = 0 \end{array} \right. \qquad [23] $$
u and v being velocity vector v components. We consider the driven cavity problem which implies the following boundary conditions:

$$ \left\{ \begin{array}{l} \underline{v}^T(x=0, y) = (0,0) \\ \underline{v}^T(x=1, y) = (0,0) \\ \underline{v}^T(x, y=0) = (0,0) \\ \underline{v}^T(x, y=1) = (x-x^2, 0) \end{array} \right. \qquad [24] $$

If all the unknown fields are approximated using standard RKPM with second degree of consistency, with a single degree of freedom by node, spurious oscillations are found in both the velocity and the pressure field. Figures 6, 7 and 8 depict both velocity components as well as the pressure field.

<figure><img src="image_6.png" /><figcaption>Figure 6. x-component of the velocity field in the driven cavity problem</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7. y-component of the velocity field in the driven cavity problem</figcaption></figure>

These results were computed by using 15x15 nodes and 900 integration points uniformly distributed.
<figure><img src="image_8.png" /><figcaption>Figure 8. Pressure field in the driven cavity problem</figcaption></figure>

To avoid this spurious effects one could try to impose different degree of reproduction for both fields: second and first degrees for the velocity and pressure respectively. However, this alternative does not remove the instabilities [3]. The use of different support sizes of the shape functions is also inefficient [3]. For this reason, we analyze the behavior associated with the use of a Hermite-RKPM approximation, being the reproduction matrix for both velocity components:

$$ \underline{R}^u(x, y) = \underline{R}^v(x, y) = \begin{pmatrix} 1 & 0 & 0 \\ x & 1 & 0 \\ y & 0 & 1 \\ x^2 & 2x & 0 \\ xy & y & x \\ y^2 & 0 & 2y \end{pmatrix} $$

being the pressure approximated using standard RKPM, whose reproduction vector results:

$$ [\underline{R}^p]^T(x, y) = (1, x, y, x^2, xy, y^2) $$

which implies 3 degrees of freedom at each node for each velocity component, and a single degree of freedom at each node for the pressure field. Figures 9, 10 and 11 depict the computed velocity and pressure fields. We can notice that the use of the H-RKPM allows stabilizing the discretization. The convergence order is only first order, possibly due to the numerical integration.
<figure><img src="image_9.png" /><figcaption>**Figure 9.** x-component of the velocity field in the driven cavity problem</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10.** y-component of the velocity field in the driven cavity problem</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 11.** Pressure field in the driven cavity problem</figcaption></figure>
**5. Conclusions**

In this paper we have proposed and applied an original Hermite reproducing kernel particle approximation (H-RKPM) for solving second order PDE. The only interest of using this kind of approximation is the possibility to reduce significantly the size of the shape functions support. Even if the support is slightly higher than the inter-nodal distance, the reproduction order can be increased by increasing the number of nodal degrees of freedom (that are the successive field derivatives).

However, its interest in mixed formulations that must verify the LBB stability condition is out of discussion. A centered formulation, as usually used in moving least squares framework, could be more suitable, simplifying the computational aspects and may be reducing the numerical errors. This is a work in progress.

**6. References**

Belytschko T., Krongauz Y., Organ D. and Fleming M., “Meshless methods: an overview and recent developments”, *Computer Methods in Applied Mechanics and Engineering*, 139, 1996, p. 3-47.

Joyot P., Trunzler J., Chinesta F., “Point collocation methods using reproducing kernel approximations for the solution of the Stokes equations”, *ECCOMAS Thematic Conference on Meshless Methods*, Lisbon, 2005.

Liu W.K., Jun S., Li S., Adee J. and Belytschko T., “Reproducing kernel particle methods for structural dynamics”, *Int. J. Numer. Methods. Eng.*, 38, 1995, p. 1655-1679.

Nayroles B., Touzot G. and Villon P., “Generalizing the finite element method: Diffuse approximation and diffuse elements”, *Journal of Computational Mechanics*, 10, 1992, p. 307-318.

Trunzler J., Joyot P. and Chinesta F., “Discontinuous derivative enrichment in RKPM”, *Lectures Notes on Computational Science and Engineering*, 43, 2005, p. 93-108.