Convergence study and optimal weight functions of an explicit particle method for the incompressible Navier--Stokes equations

TL;DR

This paper analyzes the convergence and accuracy of an explicit particle method for incompressible Navier--Stokes equations, proposing optimal weight functions and validating through numerical experiments including free surface flows.

Contribution

It introduces a convergence analysis, derives optimal weight functions for generalized approximate operators, and demonstrates improved accuracy in flow simulations.

Findings

Convergence of the explicit particle method is confirmed numerically.

Optimal weight functions reduce truncation errors.

Method effectively simulates free surface flows like dam break.

Abstract

To increase the reliability of simulations by particle methods for incompressible viscous flow problems, convergence studies and improvements of accuracy are considered for a fully explicit particle method for incompressible Navier--Stokes equations. The explicit particle method is based on a penalty problem, which converges theoretically to the incompressible Navier--Stokes equations, and is discretized in space by generalized approximate operators defined as a wider class of approximate operators than those of the smoothed particle hydrodynamics (SPH) and moving particle semi-implicit (MPS) methods. By considering an analytical derivation of the explicit particle method and truncation error estimates of the generalized approximate operators, sufficient conditions of convergence are conjectured.Under these conditions, the convergence of the explicit particle method is confirmed by numerically comparing errors between exact and approximate solutions. Moreover, by focusing on the truncation errors of the generalized approximate operators, an optimal weight function is derived by reducing the truncation errors over general particle distributions. The effectiveness of the generalized approximate operators with the optimal weight functions is confirmed using numerical results of truncation errors and driven cavity flow. As an application for flow problems with free surface effects, the explicit particle method is applied to a dam break flow.

Tables

Table 1. Table 1. Relative errors in space and time for the cases with and without the pressure recalculation ( 19 )

	(a) velocity
	without (19)	with (19)	without/with
(G-s)	$0.941$	$0.022$	$41.68$
(S-c)	$0.645$	$0.030$	$21.29$
(S-q)	$164.784$	$0.034$	$4704.57$
(S-w)	$0.686$	$0.028$	$24.07$
(M)	$0.531$	$0.034$	$15.30$

Table 2. Table 2. Convergence rates of errors obtained when Δ ​ x = 0.005 Δ 𝑥 0.005 \Delta x=0.005 and 0.0025 0.0025 0.0025

	(a) velocity
	$m=1$	$m=2$	$m=3$	$m=4$
(G-s)	-1.43	2.13	1.84	1.71
(S-c)	-1.73	2.18	1.84	1.69
(S-q)	0.37	2.14	1.84	1.64
(S-w)	-0.29	2.07	1.82	1.69
(M)	-2.18	-0.42	-0.78	0.79

Table 3. Table 3. Truncation errors of approximate Laplace operators with h = 2.1 ​ Δ ​ x , 2.6 ​ Δ ​ x , 3.1 ​ Δ ​ x ℎ 2.1 Δ 𝑥 2.6 Δ 𝑥 3.1 Δ 𝑥 h=2.1\Delta x,2.6\Delta x,3.1\Delta x and ϵ = 0.0 , 0.25 , 0.5 italic-ϵ 0.0 0.25 0.5 \epsilon=0.0,0.25,0.5 .

	$h=2.1\Delta x$
	${ϵ}_{\max }=0$	${ϵ}_{\max }=0.25$	${ϵ}_{\max }=0.5$
(G-s)	0.0532	0.6609	1.2014
(S-c)	0.0191	1.0407	1.7955
(S-q)	0.0994	1.7798	2.9894
(S-w)	0.0447	1.1538	1.9976

Table 4. Table 4. Relative errors between the reference solutions (FDM) and particle methods.

	$Re=100$			$Re=1000$
	(a)	(b)	(c)	(a)	(b)	(c)
(G-s)	0.0667	0.0484	0.0486	0.0975	0.0904	0.0877
(S-c)	0.7546	0.0513	0.0914	0.9987	0.0786	0.1068
(S-q)	1.0000	0.9206	0.0816	0.9987	0.9916	0.9998
(S-w)	0.9809	0.1206	0.0592	1.0002	0.0914	0.0835
(M)	0.1326	0.1299	0.0698	1.0757	0.9397	1.2913

Table 5. Table 5. Relative errors between experiments and particle methods at sensors.

$H_{\mathrm{dam}}$	1	2	2L	3	4
0.3	0.3829	0.2374	0.2363	0.2161	0.1834
0.6	0.2939	0.2307	0.2275	0.1953	0.1998

Equations

$$\frac{\mathrm{D}u}{\mathrm{D}t}=-\frac{1}{\rho }\nabla p+\nu \Delta u+f,$$

$$\nabla \cdot u=0,$$

$$u=u_0,$$

$$u=u_{\Gamma },$$

$${\Omega }_H:=\left\{x\in {\mathbb{R}}^d|\exists y\in \Omega \text{mbox}s.t.|x-y|<H\right\}.$$

$${\mathcal{X}}_N$$

$${\mathcal{V}}_N$$

$$\mathcal{W}:=\left\{w:[0,\infty )\to \mathbb{R}|w(r)>0(0<r<1),w(r)=0(r\ge 1)\right\}.$$

$$w_h(r):=\frac{1}{h^d}w\left(\frac{r}{h}\right),r\in [0,\infty ).$$

$$C_k(w):={\int }_{{\mathbb{R}}^d}|x{|}^{k}w(|x|)\mathrm{d}x.$$

$${\Pi }_h{\varphi }_i$$

$${\nabla }_h{\varphi }_i$$

$${\Delta }_h{\varphi }_i$$

$$C_{0,h}(w):=\frac{|{\Omega }_H|}{N}\sum _{z\in {\mathbb{Z}}^d}{w}_h(|{\Omega }_H{|}^{1/d}{N}^{-1/d}\left|z\right|),$$

$${\Lambda }^k(S):=\left\{i=1,2,\dots ,N|x_i^k\in S\right\}.$$

$${\hat{\Pi }}_h^k{\varphi }_i$$

$${\nabla }_{h,+}^k{\varphi }_i$$

$$\left\{\begin{array}[]{@{\,}r@{\,}c@{\,}ll}\dfrac{\mathrm{u}_{i}^{\ast\!\!\;,\,k+1}-\mathrm{u}_{i}^{k}}{\tau}&=&\nu\Delta_{h}^{k}\mathrm{u}_{i}^{k}+f(x_{i}^{k},t^{k}),&i\in\Lambda^{k}(\Omega),\\ \mathrm{u}_{i}^{\ast\!\!\;,\,k+1}&=&u_{\Gamma\!\!\;,\,H}(x_{i}^{k},t^{k}),&i\in\Lambda^{k}(\Gamma_{H});\end{array}\right.$$

$$\left\{\begin{array}[]{@{\,}r@{\,}c@{\,}ll}x_{i}^{\ast\!\!\;,\,k+1}&=&x_{i}^{k}+\tau\,\mathrm{u}_{i}^{\ast\!\!\;,\,k+1},&i\in\Lambda^{k}(\Omega),\\ x_{i}^{\ast\!\!\;,\,k+1}&=&x_{i}^{k},&i\in\Lambda^{k}(\Gamma_{H});\end{array}\right.$$

$${\mathrm{p}}_i^{\ast ,k+1}=\frac{\rho }{{\epsilon }^2}\left(\frac{1}{C_{0,h}(w)}\sum _{j=1}^N{\omega }_j{w}_h(|x_j^{\ast ,k+1}-x_i^{\ast ,k+1}|)-1\right),i\in {\Lambda }^k({\Omega }_H);$$

$$\left\{\begin{array}[]{@{\,}r@{\,}c@{\,}ll}x_{i}^{k+1}&=&\displaystyle x_{i}^{\ast\!\!\;,\,k+1}-\frac{\tau^{2}}{\rho}\nabla_{h,+}^{\ast\!\!\;,\,k+1}\mathrm{p}_{i}^{\ast\!\!\;,\,k+1},&i\in\Lambda^{k}(\Omega),\\ x_{i}^{k+1}&=&x_{i}^{\ast\!\!\;,\,k+1},&i\in\Lambda^{k}(\Gamma_{H}),\end{array}\right.$$

$${\mathrm{p}}_i^{k+1}={\hat{\Pi }}_h^{k+1}{\mathrm{p}}_i^{\ast ,k+1},i\in {\Lambda }^k({\Omega }_H).$$

$$\left\{\begin{array}[]{@{\,}r@{\,}c@{\,}ll}\dfrac{\mathrm{u}_{i}^{k+1}-\mathrm{u}_{i}^{\ast\!\!\;,\,k+1}}{\tau}&=&-\dfrac{1}{\rho}\nabla_{h}^{k+1}\mathrm{p}_{i}^{k+1},&i\in\Lambda^{k}(\Omega),\\ \mathrm{u}_{i}^{k+1}&=&u_{\Gamma\!\!\;,\,H}(x_{i}^{k+1},t^{k+1}),&i\in\Lambda^{k}(\Gamma_{H}).\end{array}\right.$$

$$B_r(x):=\{y\in {\mathbb{R}}^d;|x-y|<r\}.$$

$$\varphi (y)=\sum _{0\le |\alpha |\le n-1}\frac{D^{\alpha }{\varphi }_i}{\alpha !}(y-x_i{)}^{\alpha }+R_{n,i}(y;\varphi ).$$

$$R_{n,i}(y;\varphi ):=$$

$$\frac{d(n-2)!}{h^{n-2}{C}_{n-2}(w)}\frac{(y-x_i{)}^{{\beta }_{k,n-2}}}{|y-x_i{|}^{n-2}}{w}_h(|y-x_i|)$$

$$\begin{gathered} \frac{d(n-2)!}{h^{n-2}{C}_{n-2}(w)}{\int }_{{\Omega }_H}\varphi (y)\frac{(y-x_i{)}^{{\beta }_{k,n-2}}}{|y-x_i{|}^{n-2}}{w}_h(|y-x_i|)\mathrm{d}y \\ =\frac{d(n-2)!}{h^{n-2}{C}_{n-2}(w)}\sum _{0\le |\alpha |\le n-1}\frac{D^{\alpha }{\varphi }_i}{\alpha !}{\int }_{{\Omega }_H}\frac{(y-x_i{)}^{\alpha +{\beta }_{k,n-2}}}{|y-x_i{|}^{n-2}}{w}_h(|y-x_i|)\mathrm{d}y+E_{i,k,n}. \end{gathered}$$

$$E_{i,k,n}:=\frac{d(n-2)!}{h^{n-2}{C}_{n-2}(w)}{\int }_{{\Omega }_H}{R}_{n,i}(y;\varphi )\frac{(y-x_i{)}^{{\beta }_{k,n-2}}}{|y-x_i{|}^{n-2}}{w}_h(|y-x_i|)\mathrm{d}y=\mathcal{O}(h^2).$$

$$\int_{\Omega_{H}}\frac{(y-x_{i})^{\alpha+\beta_{k,n-2}}}{|y-x_{i}|^{n-2}}w_{h}(|y-x_{i}|){\rm\,d}y\\ =\begin{cases}0,\qquad&\mbox{one or more elements of $\alpha+\beta_{k,n-2}$ are odd},\\ \dfrac{h^{n-2}C_{n-2}(w)}{d},\qquad&\alpha+\beta_{k,n-2}=\beta_{k,2},\\ C_{0}(w),\qquad&n=2,\alpha=0,\end{cases}$$

$${\varphi }_i=\frac{1}{C_0(w)}{\int }_{{\Omega }_H}\varphi (y)w_h(|y-x_i|)\mathrm{d}y+\mathcal{O}(h^2)$$

$${(\nabla {\varphi }_i)}^{(k)}=\frac{d}{h{C}_1(w)}{\int }_{{\Omega }_H}\{\varphi (y)-{\varphi }_i\}\frac{{(y-x_i)}^{(k)}}{|y-x_i|}{w}_h(|y-x_i|)\mathrm{d}y+\mathcal{O}(h^2).$$

$$\sum _{k=1}^d\frac{(y-x_i{)}^{{\beta }_{k,2}}}{|y-x_i{|}^2}=1$$

$$\sum_{k=1}^{d}\int_{\Omega_{H}}(y-x_{i})^{\alpha}w_{h}(|y-x_{i}|){\rm\,d}y=\begin{cases}0,\qquad&\mbox{one or more elements of $\alpha$ are odd},\\ \dfrac{h^{2}C_{2}(w)}{d},\qquad&\mbox{$|\alpha|=2$ and all elements of $\alpha$ are even},\end{cases}$$

Full text

1. Introduction

Particle methods, such as the smoothed particle hydrodynamics (SPH) [6, 18, 25] and moving particle semi-implicit (MPS) [14, 13, 29] methods, discretize partial differential equations based on particles distributed in domains and basis functions referred to as weight functions corresponding to each particle. These particle methods do not require mesh generation; therefore, they are appropriate for problems that include large deformations or damages, e.g., collapses [20], brittle solids [3], and Navier–Stokes equations under free surface effects [14, 16, 19, 30]. In particular, explicit particle methods for Navier–Stokes equations have been widely used for large-scale problems, such as tsunami run-up [4, 22], because of their simple implementation, which can also be done using parallel computing.

Representative examples of explicit particle methods for the incompressible Navier–Stokes equations include the weekly compressible SPH (WCSPH) [19, 21] and the explicit MPS (E-MPS) [23, 28] methods. WCSPH is characterized as an explicit particle method that uses approximate differential operators of SPH for spatial discretization and evaluating pressure given an equation of state for compressible flow. In contrast, E-MPS is characterized as an explicit particle method that uses approximate differential operators of MPS for spatial discretization and evaluating pressure given the local density of a number of particles. In previous studies [19, 21, 23, 28], both these methods have been validated by comparing their numerical results with experimental results. However, in order to ensure reliability for problems on a non-experimental scale or for large-scale computations such as those involving tsunamis, numerical analyses of particle methods, such as convergence studies, are indispensable. Although there are a few mathematical analyses for particle methods or related methods [1, 2, 8, 9, 10, 26], their results do not directly apply to explicit particle methods. Therefore, we focus on convergence studies for explicit particle methods for the incompressible Navier–Stokes equations.

To describe a mathematical convergence for the incompressible Navier–Stokes equations, in this study, we configure an explicit particle method without physical parameters and assumptions in a manner similar to previous literature [19, 21, 28]. Then, we introduce a penalty problem that theoretically converges to the incompressible Navier–Stokes equations and derive the explicit particle method by discretizing the penalty problem based on mathematical theory alone. In this spatial discretization, we use generalized approximate operators, which are defined as a wider class of approximate operators for particle methods of SPH and MPS. Because of this discretization, computational procedures of the explicit particle methods closely resemble that of E-MPS. Furthermore, for the explicit particle method, we conjecture sufficient conditions of convergence based on its analytical derivation and truncation error estimates for the generalized approximate operators. Under these sufficient conditions, we confirm the convergence of the explicit particle method by computing errors between numerical solutions and exact solutions in the Taylor–Green vortex.

Moreover, to improve the accuracy of the explicit particle method, we consider an optimization of the discrete parameters based on the truncation error estimates of the generalized approximate operators [8, 11, 12]. In particular, defining the generalized approximate operators as a wider class of those used in particle methods enables us to consider an optimization of discrete parameters without imposed constraint conditions in each method. Thus, using truncation errors based on particle distributions as the objective function, we introduce an optimization problem for weight functions of the generalized approximate operators. The effects of weight functions obtained as solutions of the optimization problem are confirmed by numerical results of truncation errors and a driven cavity flow.

Furthermore, to confirm that the explicit particle method can be applied to more realistic problems, we develop it for flow problems under free surface effects. In the case of the original procedure of the explicit particle method, pressure around a free surface are evaluated as much lower than that in the inner domain of the fluid, because of the lack of particles. In addition, clustering of particles around free surface using pressure gradients causes unstable motion. Therefore, by modifying the procedure of evaluating pressure and its gradient, we ensure stable simulations of flow problems under free surface effects. Moreover, we apply the explicit particle method with these modifications to a dam break flow and compare the obtained numerical and experimental results.

2. Explicit particle method for incompressible Navier–Stokes equations

In this section, we present the formulation of the governing equations and approximate operators, which are used for spatial discretization in our study; furthermore, we introduce an explicit particle method for incompressible Navier–Stokes equations.

2.1. Governing equation

Let $\mathbb{R}$ be the set of real numbers. Let $\Omega$ be a bounded domain in $\mathbb{R}^{d}\,(d=2,3)$ with a smooth boundary $\Gamma$. We consider the incompressible Navier–Stokes equations as follows:

[TABLE]

where $u:\Omega\times(0,T)\rightarrow\mathbb{R}^{d}$, $p:\Omega\times(0,T)\rightarrow\mathbb{R}$, $\rho>0$, $\nu>0$, $f:\Omega\times(0,T)\rightarrow\mathbb{R}^{d}$, $u_{0}:\Omega\rightarrow\mathbb{R}^{d}$, and $u_{\Gamma}:\Gamma\times(0,T)\rightarrow\mathbb{R}^{d}$ denote velocity, pressure, density, kinematic viscosity, body force, initial velocity, and boundary velocity of the fluid, respectively. Furthermore, ${\rm D}/{\rm D}t$ denotes the material derivative defined as ${\rm D}/{\rm D}t\mathrel{\mathop{:}}=\partial/\partial t+u\cdot\nabla$. The unknown values are velocity $u$ and pressure $p$. We assume the uniqueness and existence of a smooth solution for the incompressible Navier–Stokes equations (1). Note that we only treat the Dirichlet boundary condition in (1) for simplicity, although we will deal with boundaries including the free surface in an applied example in Section 5.

2.2. Generalized approximate operators

We introduce approximate operators for spatial discretization of an explicit particle method. To ensure generality, we use different formulations from those used in the specific cases of SPH and MPS.

For a fixed positive number $H$ and domain $\Omega\subset\mathbb{R}^{d}$, an expanded domain $\Omega_{H}$ is defined as

[TABLE]

Let $\Gamma_{H}\mathrel{\mathop{:}}=\Omega_{H}\setminus\Omega$. Let $\mathbb{N}$ be the set of positive integers. For $N\in\mathbb{N}$, we define a particle distribution $\mathcal{X}_{N}$ and particle volume set $\mathcal{V}_{N}$ as

[TABLE]

respectively. Here, $\left|\Omega_{H}\right|$ indicates the volume of $\Omega_{H}$. We refer to $x_{i}\in\mathcal{X}_{N}$ and $\omega_{i}\in\mathcal{V}_{N}$ as a particle and particle volume, respectively. Figure 1 shows an example of a particle distribution.

We define a function space $\mathcal{W}$ as

[TABLE]

We refer to $w\in\mathcal{W}$ as a reference weight function. The influence radius $h$ is a real number satisfying $\min\{|x_{i}-x_{j}|\mid i\neq j\}<h<H$. For the reference weight function $w$ and influence radius $h$, a weight function $w_{h}:[0,\infty)\rightarrow\mathbb{R}$ is defined as

[TABLE]

We refer to the domain $\{y\in\mathbb{R}^{d};\,|y-x_{i}|<h\}$ as an influence domain for particle $x_{i}$; in addition, we refer to particles in the influence domain for particle $x_{i}$ as the neighbor particles of particle $x_{i}$. For an integer $k$ and function $w:[0,\infty)\rightarrow\mathbb{R}$, we define $C_{k}(w)$ as

[TABLE]

Set discrete parameters $(\mathcal{X}_{N},\mathcal{V}_{N},h)$ and reference weight functions $w^{\Pi},w^{\nabla},w^{\Delta}\in\mathcal{W}$. Then, for $\phi:\mathcal{X}_{N}\rightarrow\mathbb{R}$, we define the interpolant $\Pi_{h}$, approximate gradient operators $\nabla_{h}$, and approximate Laplace operator $\Delta_{h}$ as

[TABLE]

respectively. Here, $\phi_{i}\mathrel{\mathop{:}}=\phi(x_{i})$, $C_{\Pi}\mathrel{\mathop{:}}=1/C_{0}(w^{\Pi})$, $C_{\nabla}\mathrel{\mathop{:}}=d/C_{1}(w^{\nabla})$, and $C_{\Delta}\mathrel{\mathop{:}}=2d/C_{2}(w^{\Delta})$.

The derivations of these operators are presented in Section 3.1. Moreover, as discussed later in Appendix B, these operators represent a wider class of approximate operators for particle methods that those in the SPH and MPS methods. Thus, we refer to these operators as generalized approximate operators. Note that different symbols in reference weight functions $w^{\Pi},w^{\nabla},w^{\Delta}$ for each differential operator are used in order to allow us to choose them arbitrarily. In Section 4.1, we discuss the optimization of weight functions for these generalized approximate operators based on truncation error estimates.

2.3. Computational procedure of the explicit particle method

We introduce an explicit particle method for the incompressible Navier–Stokes equations. Before introducing this method, we introduce some notations used in our study. Let $u_{0\!\!\;,\,H}:\Omega_{H}\rightarrow\mathbb{R}^{d}$ and $u_{\Gamma\!\!\;,\,H}:\Gamma_{H}\times[0,T]\rightarrow\mathbb{R}^{d}$ be expanded functions of the initial and boundary velocities, respectively, that satisfy $u_{0\!\!\;,\,H}|_{\Omega}=u_{0}$, $u_{\Gamma\!\!\;,\,H}|_{\Gamma}=u_{\Gamma}$, and $u_{\Gamma\!\!\;,\,H}|_{t=0}=u_{0\!\!\;,\,H}|_{\Gamma_{H}}$. In addition, let $\tau>0$ be the time step. Further, let $K$ be the total number of time steps defined by $K\mathrel{\mathop{:}}=\lfloor T/\tau\rfloor$, where $\lfloor a\rfloor$ denotes the greatest integer that is less than or equal to $a$; this symbol is known as the Gauss symbol. For $k=0,1,\dots,K$, the $k$th time $t^{k}$ is defined as $t^{k}\mathrel{\mathop{:}}=k\,\tau$. Let $\mathcal{X}_{N}^{k}$ and $x_{i}^{k}$ be a particle distribution and an $i$th particle in that distribution at $t^{k}$, respectively. Let $\mathcal{X}_{N}^{\ast\!\!\;,\,k}$ and $x_{i}^{\ast\!\!\;,\,k}$ be a tentative particle distribution and a tentative $i$th particle in that distribution at $t^{k}$, respectively. For $w\in\mathcal{W}$, we define $C_{0,h}(w)$, which is an approximation of $C_{0}(w)$, as

[TABLE]

where $\mathbb{Z}$ is the set of integers. For $S\subset\mathbb{R}^{d}$, let $\Lambda^{k}(S)$ be an index set of particles in $S$:

[TABLE]

We denote $\Pi_{h}$, $\nabla_{h}$, and $\Delta_{h}$ by replacing $x_{i}\in\mathcal{X}_{N}$ with $x_{i}^{k}\in\mathcal{X}_{N}^{k}$ as $\Pi_{h}^{k}$, $\nabla_{h}^{k}$, and $\Delta_{h}^{k}$, respectively. For $\phi:\mathcal{X}_{N}\rightarrow\mathbb{R}$, we define a modified interpolant $\widehat{\Pi}_{h}^{k}$ and an additional approximate gradient operator $\nabla_{h,+}^{k}$ as

[TABLE]

respectively. Note that the modified interpolant $\widehat{\Pi}_{h}^{k}$ corresponds to an interpolant used for a Shepard filter [24, 27].

The computational procedure of the explicit particle method involves the following steps. Set the discrete parameters as follows: $H>0$, initial particle distribution $\mathcal{X}_{N}^{0}$, particle volume set $\mathcal{V}_{N}$, reference weight function $w$, influence radius $h\leq H/2$, parameter $\varepsilon>0$, and time step $\tau$. Set the initial approximate velocity $\mathrm{u}_{i}^{0}\,(i=1,2,\dots,N)$ as $\mathrm{u}_{i}^{0}=u_{0\!\!\;,\,H}(x_{i}^{0})$. Then, for $k=0,1,\dots,K-1$, the approximate solution $(\mathrm{u}_{i}^{k+1},\mathrm{p}_{i}^{k+1})\,(i=1,2,\dots,N)$ is solved using the following steps:

Step 1: Compute a predictor of velocity $\mathrm{u}_{i}^{\ast\!\!\;,\,k+1}$ as follows:

[TABLE]

Step 2: Compute a tentative particle position $x_{i}^{\ast\!\!\;,\,k+1}$ as follows:

[TABLE]

Step 3: Compute a tentative pressure $\mathrm{p}_{i}^{\ast\!\!\;,\,k+1}$ as follows:

[TABLE]

Step 4: Update the particle position $x_{i}^{k+1}$ as follows:

[TABLE]

where $\nabla_{h,+}^{\ast\!\!\;,\,k+1}$ is the gradient operator $\nabla_{h,+}^{k+1}$ wherein $\{x_{i}^{k+1}\}$ is replaced with $\{x_{i}^{\ast\!\!\;,\,k+1}\}$;

Step 5: Evaluate the pressure $\mathrm{p}_{i}^{k+1}$ as follows:

[TABLE]

Step 6: Evaluate the velocity $\mathrm{u}_{i}^{k+1}$ as follows:

[TABLE]

The flowchart of the explicit particle method is shown in Figure 2.

Because the pressure $\mathrm{p}_{i}^{k+1}$ is evaluated based on the density of neighbor particles, the explicit particle method is similar to E-MPS [28]. In the next section, we derive the explicit particle method and consider its errors; in addition, we show the convergence of the explicit particle method numerically.

3. Convergence study

In order to confirm the convergence of the explicit particle method, we conjecture conditions of convergence by considering the truncation error estimates of generalized approximate operators and the derivation of the explicit particle method. Moreover, we show the convergence of the explicit particle method using numerical results. Note: See Appendix A for computational rules of the multi-index and definitions of functional spaces and their norms.

3.1. Derivation of generalized approximate operators

In order to estimate truncation errors, we present the derivations of the generalized approximate operators in Section 2.2. Let $x_{i}\in\mathcal{X}_{N}\cap\Omega$. Let $B_{r}(x)$ be an open ball with the center at $x$ and radius $r$:

[TABLE]

Then, by Taylor expansion, for $y\in B_{h}(x_{i})\setminus\{x_{i}\}$ and $v\in C^{n}(\overline{\Omega}_{H})~{}(n\in\mathbb{N})$, we have

[TABLE]

Here, $\alpha$ is a multi-index and $R_{n,i}(y;\phi)$ is the residual given by

[TABLE]

For $k=1,2,\dots,d$ and nonnegative integer $l$, let $\beta_{k,l}$ be a multi-index such that the $k$th element is $l$, while the others are [math]. For $n=2,3,4$ and $k=1,2,\dots,d$, multiplying both the sides of (23) by

[TABLE]

and integrating it over $\Omega_{H}$, we get

[TABLE]

Here, $E_{i,k,n}$ is

[TABLE]

By considering that

[TABLE]

for (25) with $n=2,3$, we obtain

[TABLE]

and

[TABLE]

Moreover, when $n=4$, by

[TABLE]

and

[TABLE]

we obtain

[TABLE]

By (4), the above integration can be approximated as

[TABLE]

Therefore, by (28) and (33), and replacing $w$ with $w^{\Pi}\in\mathcal{W}$, we derive the generalized interpolant (8) as follows:

[TABLE]

By (29) and (33), and replacing $w$ with $w^{\nabla}\in\mathcal{W}$, we derive the generalized approximate gradient operator (9) as follows:

[TABLE]

Moreover, by (32) and (33), and replacing $w$ with $w^{\Delta}\in\mathcal{W}$, we derive the generalized approximate Laplace operator (10) as follows:

[TABLE]

The generalized approximate operators can be used as approximate operators of the conventional particle methods such as SPH and MPS by selecting the parameters of the generalized approximate operators appropriately; this is discussed further in Appendix B. Therefore, approximate operators of conventional particle methods can be derived using the abovementioned method.

3.2. Truncation errors of generalized approximate operators

We analyze the truncation errors of the generalized approximate operators using their derivations. Let us consider a truncation error estimate of the generalized approximate Laplace operator (10). We assume $x_{i}\in\mathcal{X}_{N}\cap\Omega$, $\phi\in C^{4}(\overline{\Omega}_{H})$, and $w^{\Delta}\in\mathcal{W}\cap C^{1}([0,\infty))$. From the derivation of the generalized approximate Laplace operator (36), we estimate its truncation error as

[TABLE]

Here,

[TABLE]

[TABLE]

Note that the estimate $\widetilde{E}_{i}=\mathcal{O}(h^{2})$ is derived from (26). Now, we estimate the error $\widehat{E}_{i}$, which consists of the integration and the numerical integration, which are the first and second terms on the right hand side of (39), respectively. For a $C^{1}$ class function $g:\Omega_{H}\rightarrow\mathbb{R}$ and generators $y_{i}\in\Omega_{H}\,(i=1,2,\dots,N)$, we assume a numerical integration for the integration of $g$ over $\Omega_{H}$ given by

[TABLE]

Here, $\sigma=\{\sigma_{i}\}_{i=1}^{N}$ is a decomposition of $\Omega_{H}$ satisfying

[TABLE]

where $\overline{\sigma}_{i}$ is the closure of $\sigma_{i}$. Then, as an estimate of the Riemann sum, we can estimate the numerical integration as

[TABLE]

Here, ${\rm rad}(\sigma_{i}):=\sup\left\{|y_{i}-z|\;\middle|\;z\in\sigma_{i}\right\}$. Furthermore, because $\sigma$ is arbitrary, we can estimate the numerical integration as

[TABLE]

From the strategy above, we introduce a decomposition of $\Omega_{H}$ as $\sigma=\{\sigma_{i}\}_{i=1}^{N}$ such that

[TABLE]

For $\sigma$, indicator $\delta_{\sigma}$ is defined as

[TABLE]

and indicator $\delta_{\infty}=\delta_{\infty}(\mathcal{X}_{N},\mathcal{V}_{N})$ is defined as

[TABLE]

Let any $\sigma=\{\sigma_{i}\}_{i=1}^{N}$ such that (44). Furthermore, we assume $\delta_{\infty}\leq h$. Then, by Taylor’s theorem, we can estimate the following:

[TABLE]

Because $\sigma$ is arbitrary, we obtain

[TABLE]

Hence, $\delta_{\infty}\leq h$ yields

[TABLE]

Consequently, by (37) and (38), and (49), we establish

[TABLE]

Let $r_{\min}$ be $r_{\min}\mathrel{\mathop{:}}=\min\{|x_{j}-x_{i}|;\,i,j=1,2,\dots,N,i\neq j\}$. If $\delta_{\infty}=\mathcal{O}(r_{\min})$, we refer to the particle distribution and particle volume as regular. Based on the definition of $\delta_{\infty}$, in the absence of extremely unfavorable conditions, such as high density particle distributions or high variance particle volumes, the particle distribution and particle volume become regular. The indicator $\delta_{\infty}$ satisfies $\delta_{\infty}=\mathcal{O}(N^{-1/d})=\mathcal{O}(r_{\min})$. By assuming the regularity of the particle distribution and particle volume, we estimate the truncation error of the generalized approximate Laplace operator as

[TABLE]

A more precise theorem to estimate truncation error has been reported in the literature[8, 11, 11].

3.3. Derivation of the explicit particle method

The explicit particle method is based on the following penalty problem for the incompressible Navier–Stokes equations:

[TABLE]

Here, $\varepsilon$ and $p_{0}$ are a penalty term in $\mathbb{R}$ and the initial pressure, respectively. Furthermore, ${\rm D}_{t,\varepsilon}\,$ denotes a material derivative defined as ${\rm D}_{t,\varepsilon}\,\mathrel{\mathop{:}}=\partial/\partial t+u_{\varepsilon}\cdot\nabla$. The unknown values include $u_{\varepsilon}:\Omega_{H}\times[0,T)\rightarrow\mathbb{R}^{d}$ and $p_{\varepsilon}:\Omega\times[0,T)\rightarrow\mathbb{R}$. (52a) is the moment equation, which is the same as (1a). Further, (52b) is based on the continuity equation for the compressible flow. If $\varepsilon=0$, we find that (52b) is equivalent to (1b). Therefore, the solution $(u_{\varepsilon},p_{\varepsilon})$ in the penalty problem (52) coincides with the solution $(u,p)$ in the original incompressible Navier–Stokes equations (1) if $\varepsilon=0$ formally. In particular, in the cases of two-dimensional spaces and partially- or full-periodic boundary conditions, the convergence of the penalty problem (52) has orders of velocity and pressure as $\mathcal{O}(\varepsilon^{2})$ and $\mathcal{O}(\varepsilon)$, respectively, which has been proved in Kreiss et al.[15].

We arbitrarily set $k=0,1,\dots,K-1$. Before deriving the discretized schemes, we define a function $\widetilde{p}^{k}_{\varepsilon}:\Omega\times[t^{k},t^{k+1})\rightarrow\mathbb{R}$ as

[TABLE]

where $w\in\mathcal{W}$ and $X^{k}_{\varepsilon}$ is the solution of the following differential equation:

[TABLE]

Then, under the assumption that $\|u_{\varepsilon}\|_{C^{1}([0,T];C^{3}(\Omega_{H}))}<\infty$, we have

[TABLE]

The proof for which is presented in Appendix C. By comparing (52b) and (55), the function $\widetilde{p}^{k}_{\varepsilon}$ yields an approximation of the pressure $p_{\varepsilon}$ at $t\in[t^{k},t^{k+1})$.

Next, we introduce a time-discretized scheme for the penalty problem (52). For $k=0,1,\dots,K$, let $(u_{\varepsilon,\tau}^{k},p_{\varepsilon,\tau}^{k})$ be a solution of this scheme at $t=t^{k}$. We set $k=1,2,\dots,K-1$ and $x\in\Omega$ in an arbitrary manner. For $y\in\Omega$, we introduce $X^{k+1}_{\varepsilon,\tau}(y)\,(\approx X^{k}_{\varepsilon}(y,t^{k+1}))$ as

[TABLE]

Because the material derivative is estimated by

[TABLE]

we discretize (52a) as

[TABLE]

To evaluate the approximate pressure $p_{\varepsilon,\tau}^{k+1}$, we introduce a tentative velocity $u_{\varepsilon,\tau}^{\ast\!\!\;,\,k+1}$ and tentative position $X^{\ast\!\!\;,\,k+1}_{\varepsilon,\tau}$ for $x\in\Omega$ as

[TABLE]

and

[TABLE]

respectively. Using these equations and (53), the approximate pressure $p_{\varepsilon,\tau}^{k+1}$ is obtained as follows:

[TABLE]

Then, by discretizing the time-discretized scheme in space, we derive the explicit particle method. Let $i=1,2,\dots,N$ such that $x_{i}^{k}\in\Omega$. First, we discretize (59) and (60) as

[TABLE]

and

[TABLE]

respectively. Using these, we discretize (61) as

[TABLE]

Then, the particle position is updated as follows:

[TABLE]

In this case, to avoid non-uniform particle distributions as discussed in Price [25], we use $\nabla_{h,+}^{\ast\!\!\;,\,k+1}$ as the gradient operator in (65). Furthermore, to eliminate noise corresponding to the density of particles, we modify the pressure calculation as follows:

[TABLE]

Then, by using pressure $\mathrm{p}_{i}^{k+1}$, we discretize (58) as

[TABLE]

Finally, using (62) and (67), we derive (20). The pressure recalculation (66) is essential to obtain stable and accurate results as shown in numerical experiments in Section 3.5.

3.4. Sufficient conditions of convergence

We conjecture the sufficient conditions of convergence for the explicit particle method by considering the deviations and truncation error estimates that were calculated in previous sections. In particular, as sufficient conditions of convergence, we require $h\rightarrow 0$ and $r_{\min}h^{-2}\rightarrow 0$ from the truncation error estimates (51) calculated in Section 3.2. In addition, because the convergence orders between the solution of the incompressible Navier–Stokes equation (1) and that of the penalty problem (52) are $\mathcal{O}(\varepsilon)$, we require $\varepsilon\rightarrow 0$. Moreover, we require $\tau h^{-1}\rightarrow 0$ and $\tau\rightarrow 0$ from the order estimates (55) and (57) obtained in Section 3.3. Thus, the summary of the above conditions is as follows: $h\rightarrow 0$, $r_{\min}h^{-2}\rightarrow 0$, $\varepsilon\rightarrow 0$, $\tau h^{-1}\rightarrow 0$, and $\tau\rightarrow 0$. In particular, when the time step $\tau$ satisfies

[TABLE]

where $\|\cdot\|_{L^{\infty}([0,T];L^{\infty}(\Omega))}$ denotes the infinity norm in space-time $\|\phi\|_{L^{\infty}([0,T];L^{\infty}(\Omega))}\mathrel{\mathop{:}}=\operatorname{ess\,\sup}\{|\phi(x,t)|;and\,x\in\Omega,t\in(0,T)\}$; the conditions for convergence are

[TABLE]

Condition (68) is based on the von Neumann stability analysis and corresponds to that suggested in Morris et al. [21] by replacing the sound speed with $\varepsilon^{-1}$. Under the conjectured sufficient conditions, we confirm the convergence of the explicit particle method numerically; this is shown in the next subsection.

3.5. Numerical convergence

We confirm the conjectured sufficient condition of convergence using the numerical results for the Taylor–Green vortex. The Taylor–Green vortex is one of the solutions of the two-dimensional incompressible Navier–Stokes equations (1) in the absence of body force ($f\equiv 0$). Let $\Omega=(0,L)\times(0,L)$. The solutions of the Taylor–Green vortex $(u=(u_{1},u_{2})^{T},p)$ are given by

[TABLE]

Here, $U$ is the velocity scale, and $Re$ is the Reynolds number defined as $Re\mathrel{\mathop{:}}=UL/\nu$. Hereafter, we set $T=0.1$, $\rho=1$, $U=1$, $L=1$, and $\nu=10^{-1}$, namely, $Re=10$. By comparing the exact solution and a numerical solution of the Taylor–Green vortex, we investigate the validity of the accuracy of the pressure recalculation (19) and convergences of the explicit particle method. It should be noted that we do not treat a comparison of accuracy for approximate operators here because the Taylor–Green vortex represents an isotopic flow and disturbances in particle distributions rarely appear in the case when the explicit particle method is used.

Before performing the numerical experiments for convergence, we confirm the computational stability and accuracy of the explicit particle method. Because the Taylor–Green vortex is periodic in space, we consider a periodic domain. In particular, we consider the following coordinate system: (${x}^{(1)},{x}^{(2)}$); ${x}^{(k)}\leftarrow{x}^{(k)}+1$ if ${x}^{(k)}\leq 0$ and ${x}^{(k)}\leftarrow{x}^{(k)}-1$ if ${x}^{(k)}\geq 1$ for $k=1,2$. Then, the particles near the boundary refer to the influence domain corresponding to the periodic boundary conditions as shown in Figure 3. Moreover, if a particle crosses over the boundary, we let the particle move according to the treatment shown in Figure 4.

Because the boundary condition is not required for the system, we do not set the parameter $H$ and expanded boundary condition $u_{\Gamma\!\!\;,\,H}$ for it. The initial particle distribution $\mathcal{X}_{N}^{0}$ is set as the square lattice with spacing $\Delta x=0.04$:

[TABLE]

Then, the number of particles is $N=(\lfloor 1/\Delta x\rfloor)^{2}=25^{2}$. Furthermore, the particle volume set $\mathcal{V}_{N}=\{\omega_{i}\}_{i=1}^{N}$ is set as

[TABLE]

We consider five sets of reference weight functions $(w^{\Pi},w^{\nabla},w^{\Delta})$;

1. (G-s)

$w^{\Pi}$, $w^{\nabla}$, and $w^{\Delta}$ are set as

[TABLE]

where $w^{\rm spike}$ is the spike function given by

[TABLE] 2. (S-c)

$w^{\Pi}$, $w^{\nabla}$, and $w^{\Delta}$ are set as

[TABLE]

where $\dot{w}^{\rm SPH}$ is the first derivative of $w^{\rm SPH}$, in which $w^{\rm SPH}$ uses the cubic B-spline defined as

[TABLE] 3. (S-q)

$w^{\Pi}$, $w^{\nabla}$, and $w^{\Delta}$ are set by (77) in which $w^{\rm SPH}$ uses the quintic B-spline defined as

[TABLE] 4. (S-w)

$w^{\Pi}$, $w^{\nabla}$, and $w^{\Delta}$ are set by (77) in which $w^{\rm SPH}$ uses the quintic Wendland function (a fifth positive definite function) defined as

[TABLE] 5. (M)

$w^{\Pi}$, $w^{\nabla}$, and $w^{\Delta}$ are set as

[TABLE]

respectively, where $w^{\rm MPS}$ is the reference weight function of MPS defined as

[TABLE]

Here, $a_{d}^{\rm cubic}$, $a_{d}^{\rm quintic}$, and $a_{d}^{\rm Wendland}$ are constants that satisfy the unity condition:

[TABLE]

As shown in Appendix B, the cases (S-c), (S-q), and (S-w) correspond to the use of approximate operators in SPH. Further, case (M) corresponds to the use of approximate operators in MPS. The influence radius is set as $h=0.124\,(=3.1\Delta x)$. In addition, we set $\varepsilon=0.1$ and $\tau=\tau_{\max}$.

Under the computational settings above, for the explicit particle method and that without the pressure recalculation (19), we compute the relative errors in space as:

[TABLE]

and the relative errors in space and time as:

[TABLE]

Here, the norms are defined as

[TABLE]

[TABLE]

Moreover, $\overline{\mathrm{p}}$ is defined by

[TABLE]

Then, $\overline{\mathrm{p}}$ satisfies the following condition:

[TABLE]

This condition corresponds to the integration condition of pressure:

[TABLE]

Figure LABEL:fig:TGvortex_error_time_history shows time histories of the relative errors; in the figure, the vertical axes are plotted on the logarithmic scale. Table 1 lists the relative errors in space and time of velocity as well as pressure. In all the cases except for (M), the errors of pressure, which first oscillate, are considerably improved with the pressure recalculation compared with the case without the pressure recalculation. Moreover, the accuracy of velocity is enhanced by improving the accuracy of pressure. In the case of (M), the accuracy of pressure is not remarkably different between the cases with and without pressure recalculation; nevertheless, the accuracy of velocity is improved with pressure recalculation, which is clear from Table 1. Consequently, we use the method involving the pressure recalculation.

Next, we investigate the convergence of approximate solutions for the explicit particle method. We set the initial particle distributions using (73) with $\Delta x=0.04,0.02,0.01,0.005,0.0025$. The particle volume set and reference weight functions are set as the same in the previous cases. For $m=1,2,3,4$, the influence radius $h$ is set as $h=C_{m}\Delta x^{1/m}$. Here, $C_{m}$ can be obtained by $C_{m}=3.1\times(0.04)^{1-1/m}$, which satisfies $h=3.1\times 0.04$ when $\Delta x=0.04$ for all $m$. We set $\tau=\tau_{\max}$ and $\varepsilon=2.5\Delta x$. Then, by assuming that particles maintain distances proportional to $\Delta x$, i.e., $r_{\min}=\mathcal{O}(\Delta x)$, we have $\varepsilon=\mathcal{O}(\Delta x)$, $h=\mathcal{O}(\Delta x^{1/m})$, and $r_{\min}h^{-2}=\mathcal{O}(\Delta x^{1-2/m})$. Therefore, the conditions (69) are satisfied when $\Delta x\rightarrow 0$ in the case that $m>2$. It should be noted that although $m=1$ is often used in practical computing because the average number of particles in the influence domain increases exponentially when $m>1$, as shown in Figure 6, we must set $m>1$ to conduct simulations with convergence.

Under these conditions, we compute the relative errors in space and time. Figure 7 shows the double logarithmic graph of the relative errors versus the influence radius $h$. Here, the slopes of the hypotenuse of the triangle in Figure 7 ($m\geq 2$) are $\mathcal{O}(h^{2})$ for (a) velocity and $\mathcal{O}(h^{(m-1)/2})$ for (b) pressure. Table 2 lists the convergence rates of velocity and pressure obtained for $\Delta x=0.005$ and $0.0025$. From Figure 7 and Table 2, we can confirm that the convergence orders of velocity and pressure with respect to the influence radius $h$ are of the second order and $(m-1)/2$th order for $m\geq 2$, respectively, except in case (M). This is because the approximate solution did not converge in the case of the approximate operators of MPS; this might be attributed to the fact that sufficient conditions of convergence were derived under assumptions of sufficiently smooth and bounded weight functions for the truncation error estimates.

4. Approaches for reducing truncation errors of the generalized approximate operators

In order to conduct more accurate simulations, we improve the accuracy of the generalized approximate operators by considering an optimization problem for weight functions derived based on their truncation error estimates. Moreover, the accurate results of the explicit particle methods (i.e., with an optimal weight function included) are confirmed by numerical truncation errors and numerical errors of cavity flow.

4.1. Optimization problem for weight functions

We derive an optimization problem of truncation errors for weight functions and its solutions. As discussed in Section 3.2, the truncation error of the generalized approximate Laplace operator is estimated by

[TABLE]

where

[TABLE]

[TABLE]

We can estimate $|\widetilde{E}_{i}|$ using $\mathcal{O}(h^{2})$ independent of particle distributions. Thus, we can estimate that $|\widehat{E}_{i}|$ represents an error based on disturbances of the particle distribution. In practical computing, it is rare for the particle distribution to become sufficiently uniform in each time step; hence, we aim to reduce the error $|\widehat{E}_{i}|$. In Section 3.2, we estimated $|\widehat{E}_{i}|$ as

[TABLE]

under the condition $w\in\mathcal{W}\cap C^{1}([0,\infty))$. Therefore, by using this term with respect to the weight function in $|\widehat{E}_{i}|$ as the objective function, we define the following optimization problem for the weight functions:

[TABLE]

In order to reduce the computational complexity, we consider solving the optimization problem within a range where the reference weight function transforms into a polynomial function. We give the reference weight function as the $n$th polynomial function:

[TABLE]

Because the condition $w\in\mathcal{W}\cap C^{1}([0,\infty))$ yields

[TABLE]

we have the conditions of the coefficients in (98):

[TABLE]

Therefore, in the case of a quadratic polynomial $n=2$, the solution of (97) is the spike function (76). When $n\geq 3$, we consider that the additional condition $\mbox{minimizes}~{}F(w)$ because we calculate $F(w)$ as

[TABLE]

By solving the minimization problem of $F(w)$ for the coefficients of polynomial functions under the constraint conditions (100), we can obtain optimal weight functions for $n\geq 3$. However, because the optimal weight functions with $n\geq 3$ depend on the spatial dimension, we use the quadratic spike function (76) to avoid such spatial dimension dependency in the subsequent numerical experiments.

4.2. Numerical results of truncation errors

In order to verify the analytical discussions presented in the previous section, we compute the numerical truncation errors of approximate Laplace operators when the disturbances of the particle distribution are changed. Furthermore, the test function is set as $v({x}^{(1)},{x}^{(2)})=\sin(2\pi({x}^{(1)}+{x}^{(2)}))$. The domain is set as $\Omega=(0,1)\times(0,1)$. Let $H=3\times 2^{-4}$. Then, the particle distribution is set as

[TABLE]

Here, $\Delta x=2^{-4}$ and $\epsilon_{ij}^{(k)}\,(k=1,2)$ is a random number satisfying $|\epsilon_{ij}|\leq\epsilon_{\max}\,(0\leq\epsilon_{\max}<1)$. Figure 8 shows examples of the particle distributions with perturbation $\epsilon_{\max}=0,0.25,0.5$ in $\Omega$.

It should be noted that this particle distribution becomes a square lattice if $\epsilon_{\max}=0$. The particle volume set $\mathcal{V}_{N}=\{\omega_{i}\}_{i=1}^{N}$ is determined by (74). The influence radius is set as $h=2.1\Delta x,2.6\Delta x,3.1\Delta x$. We consider the following four reference weight functions:

1. (G-s)

$w^{\Delta}$ is set as the quadratic spike function (76); 2. (S-c)

$w^{\Delta}$ is set as

[TABLE]

where $w^{\rm SPH}$ is the cubic B-spline (78). 3. (S-q)

$w^{\Delta}$ is set as (103) where $w^{\rm SPH}$ is the quintic B-spline (79); 4. (S-w)

$w^{\Delta}$ is set as (103) where $w^{\rm SPH}$ is the quintic Wendland function (80).

Figure 9 shows the graphs for the relative truncation error

[TABLE]

versus the perturbation $\epsilon_{\max}$ when $h=2.1\Delta x,2.6\Delta x,3.1\Delta x$. Table 3 lists the relative truncation errors with $h=2.1\Delta x,2.6\Delta x,3.1\Delta x$ and $\epsilon_{\max}=0.0,0.25,0.5$. From Figure 9 and Table 3, we can confirm that the truncation errors increase as perturbation $\epsilon_{\max}$ increases and influence radius rate $h$ decreases. In all the cases, though the truncation error of the generalized approximate Laplace operator with the spike function is larger than that of the conventional Laplace operators for uniform particle distributions ($\epsilon_{\max}=0$), the truncation error becomes smaller for general particle distributions ($\epsilon_{\max}>0$). Therefore, we confirmed that truncation errors can be effectively reduced for general particle distributions using the generalized Laplace operator with the spike function. Later, in Section 4.3, we confirm whether the generalized approximate operators with the spike function are also valid for fluid simulations.

4.3. Numerical results of driven cavity flow

In order to investigate whether the generalized approximate operators with the optimal weight function are also effective for a flow problem, we apply the explicit particle method to a driven cavity and compare errors in the cases when five pairs of weight functions are used. The driven cavity flow is a viscous flow problem in a rectangular domain with Dirichlet boundary conditions. One side of the boundaries flows in a tangential direction, while the other sides are wall boundaries. In the case of a square domain, the driven cavity flow can be denoted solely on the basis of the Reynolds number $Re=LU/\nu$, where $L$ and $U$ are the length of one side of the domain and velocity on the driven boundary, respectively. Hereafter, we consider $Re=100$ and $1000$.

We consider the driven cavity flow in the square domain $\Omega=(0,1)\times(0,1)$. We denote the velocity as $u=(u_{1},u_{2})^{T}$. The initial conditions are given by

[TABLE]

while the boundary conditions are given by

[TABLE]

Furthermore, zero gravity is assumed ($f=0$).

We set the parameters as follows. Set $H=0.1$, $u_{0\!\!\;,\,H}=0$ in $\Omega_{H}$ and $u_{\Gamma\!\!\;,\,H}=(1,0)\mbox{~{}in~{}}\{(x,t)\in\Omega_{H}\times(0,T);\,{x}^{(2)}\geq 1\}$, $=(0,0)\mbox{~{}in~{}}\{(x,t)\in\Omega_{H}\times(0,T);\,{x}^{(2)}<1\}$ . The initial particle distribution $\mathcal{X}_{N}^{0}$ is set as a square lattice with spacing $\Delta x$. Here, $\Delta x$ is set as $\Delta x=0.005$ and $\Delta x=0.0025$ when $Re=100$ and $Re=1000$, respectively. Note that the particles are distributed in $\Omega_{H}=(-H,1+H)\times(-H,1+H)$, and the particle distributions outside of the wall boundary correspond to well-known dummy particles [30]. We consider the same five pairs of approximate operators used in Section 3.5. We consider the same three cases of influence radii as in the numerical experiments in Section 4.2: (a) $h=2.1\Delta x$; (b) $h=2.6\Delta x$; (c) $h=3.1\Delta x$. We set $\varepsilon=0.1$ and $\tau=\tau_{\max}$. Under these conditions, we compute the two-dimensional driven cavity flow and compare velocity profiles in the vertical direction on the lines $x=0.5$ with the reference solutions, which are the numerical results of the higher-order finite difference method by Ghia et al. [5].

Figure 10 shows the velocity profiles of the two-dimensional driven cavity flow at $Re=100$ and $Re=1000$. The boxes in Figure 10 show the vertical velocity $u_{2}^{\rm\,FDM}(x^{\rm\,FDM}_{j})$ of the observation point $x^{\rm\,FDM}_{j}$ in the results of Ghia et al. [5]. Table 4 lists the errors of velocity measured using the following discrete $L^{2}$ norm in space:

[TABLE]

where $\Delta x_{j}^{\rm\,FDM}=x_{j}^{\rm\,FDM}-x_{j-1}^{\rm\,FDM}$, and

[TABLE]

From the Figure 10 and Table 4, it is clear that the velocity becomes stable as the influence radius increases. In particular, because case (G-S) has a solution even in the case of (a) $h=2.1\Delta x$, the explicit particle method using the generalized approximate operators with the spike weight function (76) is more robust to the influence radius than that with other weight functions. This result is consistent with that of the generalized approximate operators with the spike weight function (76) being more accurate for non-uniform particle distributions than other operators for the truncation error estimates, as discussed in Sections 4.1–4.2. Therefore, we confirm that the generalized approximate operators with the spike function (76) are also effective for a flow problem.

5. Application for incompressible viscous flow problems under free surface effects

In order to confirm that the explicit particle method is applicable to realistic problems, we develop the explicit particle method for flow problems under free surface effects. We introduce modifications for pressure evaluation and pressure gradient to avoid clustering of particles in and around a free surface. Moreover, we apply the modified explicit particle method to a dam break flow and compare the numerical results with experimental results.

5.1. Treatment of free surface

Because particles around the free surface do not have a sufficient number of particles in their influence domain, approximate operators on these particles do not behave appropriately. In particular, particles around the free surface come close or collide with each other in the case of the original scheme. For this reason, the tentative densities on particles around the free surface are evaluated to be considerably lower than that an inner particle. Consequently, the tentative pressure on these particles become negative per (17); then, retraction forces are experienced owing to the pressure gradient in (20). In order to solve this problem, we have to modify evaluations of the tentative density and tentative pressure. Thus, we modify (17) and (20). In order to avoid obtaining negative pressures, we modify (17) as

[TABLE]

Moreover, when the original pressure gradient is used in (20), a non-physical force develops in the tangential direction of the free surface because of the lack of a sufficient number of particles in and around the free surface. Therefore, we modify (20) as

[TABLE]

Only with the modifications above, we attain a stable and accurate simulation of a dam break flow in the next section; however, we observe strange particle motions around the free surface. Therefore, we add the collision methods [29] used in E-MPS, which modify the particle distributions to maintain their momentum, and this was confirmed to solve the problem.

5.2. Dam break flow

The dam break flow is a flow problem in which a water column on one side of a tank collapses because of gravity. Because a considerable amount of experimental data, including flow tip speeds, wave height history, and wall pressure distributions, have been collected in the literature [17, 7], changes in free surface geometry and pressure distributions in the numerical results can be confirmed.

We consider the hydraulic experiment by Lobovskỳ et al. [17] as shown in Figure 11. In this experiment, five pressure sensors are set on the opposite side of the water column. As shown in the left part of Figure 11, the five pressure sensors labeled as 1, 2, 2L, 3, and 4. In particular, their coordinates are $(0,0.075,0.003)$, $(0,0.075,0.015)$, $(0,0.0375,0.015)$, $(0,0.075,0.03)$, and $(0,0.075,0.08)$, respectively, from the origin $o$. The height of the water column $H_{\rm dam}$ is $0.3$ m or $0.6$ m.

We set the end time $T$ as $T=1.3$. Furthermore, we set the remaining parameters as follows. The initial particle distribution in the flow domain is set as a cubic lattice with spacing $\Delta x=3.0\times 10^{-3}$ m in the water column. Moreover, we set particles on a cubic lattice on the outer domain whose distances from the wall are less than $H=5.2\Delta x$. Note that the particle distributions outside of the wall boundary correspond to well-known dummy particles [30]. The velocity of the particles outside the domain are set as zero. Then, we set $\varepsilon=0.05$, $h=2.6\Delta x$, and $\tau=\tau_{\max}$. Under these conditions, we compute the dam break flow and compare the pressure at the sensors. Here, the pressures at the sensors are computed using the pressure of the nearest particle on the wall boundary from these sensors, i.e., the numerical pressure of Sensor $l$ at $t=t^{k}$ is computed as

[TABLE]

where $X_{l}$ is the position of Sensor $l$.

Figure 12 shows the pressure distributions of the explicit particle method when $H_{\rm dam}=0.3$. Figure 13 shows pressure histories of the experimental and numerical results at each sensor when $H_{\rm dam}=0.3$ m and $0.6$ m. Table 5 lists relative errors of pressure in a discrete $L^{2}$ norm in time as

[TABLE]

for Sensor $l$. Here, $P^{\rm\,ex}_{l}(t^{k})$ is the observed pressure for Sensor $l$ at $t=t^{k}$. From Figure 12, we can observe smooth pressure distributions. Moreover, from Figure 13 and Table 5, we can obtain the numerical results based on the experiment results. These numerical results show that the explicit particle method is applicable for flow problems under free surface effects.

6. Conclusion

We conducted a convergence study for an explicit particle method for the incompressible Navier–Stokes equations. The explicit particle method is based on a penalty problem of the incompressible Navier–Stokes equations, which was derived using the mathematical discretization procedure. Moreover, the explicit particle method uses generalized approximate operators, which was introduced as a wider class of approximate operators than those used in SPH and MPS for spatial discretization. By referring to the convergence orders of the penalty problem and orders of the residual appearing in the derivation process as well as truncation errors of the generalized approximate operators, we conjectured sufficient conditions of convergence for the explicit particle method. The convergence with these sufficient conditions was confirmed using numerical results of the Taylor–Green vortex; in particular, these numerical convergence orders of velocity and pressure with respect to the influence radius $h$ were $\mathcal{O}(h^{2})$ and $\mathcal{O}(h^{(m-1)/2})$ with $m\geq 2$, respectively, where $m$ is a parameter determining the ratio of increase of neighbor particles in influence.

Next, we optimized the reference weight functions considering the decreasing truncation errors of the generalized approximate operators for non-uniform particle distributions. Because the generalized approximate operators were defined as the generalization of those in conventional particle methods, we could set an optimization problem under wider conditions of parameters than those imposed in conventional particle methods. Consequently, the reference weight functions that served as the solution to the optimization problem were different from reference weight functions typically used in conventional particle methods; improvements of accuracy for non-uniform particle distributions were observed through numerical results of the truncation errors and driven cavity flow.

Finally, we developed the explicit particle method for incompressible Navier–Stokes equations with free surface effects. We modified the evaluation of pressure and approximate gradient operator in the explicit particle method to prevent the particle concentrations around the free surface becoming dense. We applied the explicit particle method with these modifications to the dam break flow and confirmed a smooth pressure distribution as well as agreement of the time histories of pressure with the experimental results.

As future work, we will investigate the stability and convergence of the particle methods mathematically. Moreover, we will develop particle methods with convergence under more practical conditions such as a that involving fixing the number of neighbor particles $(m=1)$.

Acknowledgment

This study was partly supported by priority project 3 for the Post-K Computer entitled “Sophisticated numerical analysis of diverse earthquake and tsunami disaster scenarios”.

Appendix A Notation

First, we summarize the computational rules of the multi-index. Let $\alpha=(\alpha_{1},\alpha_{2},\dots,\alpha_{d})$ be the $d$th multi-index. For a vector $x\in\mathbb{R}^{d}$, we denote the $k$th element of $x$ as ${x}^{(k)}$. Then, that operations for the multi-index are defined by

[TABLE]

Let $D^{\alpha}$ be the differential operator defined by

[TABLE]

where $D^{\alpha}v=v$ if $|\alpha|=0$.

Next, we introduce some functional spaces. For a set $S\subset\mathbb{R}^{d}\,(d\in\mathbb{N})$, let $C(\overline{S})$ be the space of real continuous functions defined in $\overline{S}$, where $\overline{S}$ is the closure of $S$. The norm of $C(\overline{S})$ is defined by

[TABLE]

For an open set $S$ and positive integer $k$, let $C^{k}(\overline{S})$ be the space of functions in $C(\overline{S})$ with derivatives up to the $k$th order. The norm of $C^{k}(\overline{S})$ is defined as

[TABLE]

Here, $\alpha$ is the multi-index. For a functional space $X(\overline{S})$, let $C([0,T];X(\overline{S}))$ be the space of functions on $\overline{S}\times[0,T]$ satisfying

[TABLE]

Appendix B Description of approximate operators in SPH and MPS using generalized approximate operators

We show that the generalized approximate operators (8)–(10) denote approximate operators in SPH and MPS if their parameters are selected appropriately. Let $w^{\rm SPH}\in\mathcal{W}$ be a reference weight function such that

[TABLE]

where $\dot{w}^{\rm SPH}$ is the first derivative of $w^{\rm SPH}$. Then, in SPH, the interpolant $\Pi_{h}^{\rm SPH}$, approximate gradient operator $\nabla_{h}^{\rm SPH}$, and approximate Laplace operator $\Delta_{h}^{\rm SPH}$ are defined as

[TABLE]

respectively. Here, $m_{j}$ and $\rho_{j}$ are positive parameters referred to as the particle mass and particle density, respectively. The particle volume set $\mathcal{V}_{N}$ is given by $\mathcal{V}_{N}=\{\omega_{i}=m_{i}/\rho_{i}\mid i=1,\dots,N\}$. Then, from (120), the generalized interpolant (8) with $w^{\Pi}=w^{\rm SPH}$ is equivalent to the interpolant of SPH (122). From

[TABLE]

the generalized approximate gradient operator (9) with $w^{\nabla}=-\dot{w}^{\rm SPH}$ is equivalent to the approximate gradient operator of SPH (123). Moreover, from (125), the generalized approximate Laplace operator (10) with

[TABLE]

is equivalent to the approximate Laplace operator of SPH (124).

Let $w^{\rm MPS}\in\mathcal{W}$ be a reference weight function defined by (82). A weight function $w^{\rm MPS}_{h}$ is set by (6). Then, in MPS, the approximate gradient operator $\nabla_{h}^{\rm MPS}$ and approximate Laplace operator $\Delta_{h}^{\rm MPS}$ are defined as

[TABLE]

respectively. Here, $n_{0}$ and $\lambda_{0}$ are parameters that depend on both $w^{\rm MPS}$ and $h$. In general, $\lambda_{0}$ is given by $\lambda_{0}=C_{2}(w^{\rm MPS}_{h})$. Then, the particle volume set $\mathcal{V}_{N}$ is given by $\mathcal{V}_{N}=\{\omega_{i}=C_{0}(w^{\rm MPS}_{h})/n_{0}\mid i=1,\dots,N\}$. Further, the generalized approximate gradient operator (9) with

[TABLE]

is equivalent to the approximate gradient operator of MPS (127). Furthermore, the generalized approximate Laplace operator (10) with $w^{\Delta}=w^{\rm MPS}$ is equivalent to the approximate Laplace operator of MPS (128) with $\lambda_{0}=C_{2}(w^{\rm MPS}_{h})$.

Appendix C Order estimates of approximate pressure

Here, we derive the order estimate (55). We assume $\|u_{\varepsilon}\|_{C^{1}([0,T];C^{3}(\Omega_{H}))}<\infty$. We arbitrarily set $k=1,2,\dots,K$. Let $t\in[t^{k-1},t^{k+1}]$. Then, by the chain rule, we have

[TABLE]

Further, by Taylor expansion and using our assumptions, we get

[TABLE]

Using the multi-indices $\alpha$ and $\beta$, we have

[TABLE]

If $\alpha=\beta$, then, by the Gauss–Green theorem and considering

[TABLE]

we have

[TABLE]

If $\alpha\neq\beta$, then, by the symmetry of the integrated function, we have

[TABLE]

Thus, we obtain

[TABLE]

Moreover, by the symmetry of the integrated function, the second term on right-hand side in (131) becomes

[TABLE]

Therefore, by (130), (131), (136), and (137), we obtain

[TABLE]