# Verification of the Born rule via direct measurement of superposition
wavefunction

Meng-Jun Hu<EMAIL_ADDRESS>Beijing Academy of Quantum Information Sciences,
Beijing 100193, China

###### Abstract

The Born rule, which is one of foundational axioms of quantum theory, states
that the probability of obtain outcome $a$ for the quantum state
$|\psi\rangle$ is determined by $P(a)=|\langle a|\psi\rangle|^{2}$. Despite
its great success in predicting the experimental outcomes, there still lacks a
direct way to verify the Born rule. Here, we show that the weak value based
direct measurement of superposition wavefunction is feasible, which can be
used to verify the Born rule directly. The plausible experiment is suggested.

##  I Introduction

For a quantum system described by the state $|\psi\rangle$, if we perform the
measurement of observable $\hat{A}$, the Born rule tells us that the
probability of obtain outcome $a$ is given by $P(a)=|\langle
a|\psi\rangle|^{2}$. The Born rule has already withstood numerical
experimental tests and there seems no question for its correctness. The
interesting and important question now is whether or not the Born rule can be
directly verified. This is possible if we can directly measure the probability
amplitude $\langle a|\psi\rangle$. For the observable to be position
$\hat{X}$, it implies the direct measurement of wavefunction $\psi(x)$. Since
the quantum state can be superposed, the best way to verify the Born rule is
to realize direct measurement of superposition wavefunction
$\psi(x_{1})+\psi(x_{2})$ and $\psi(x_{1}),\psi(x_{2})$ respectively, and then
check if they satisfy the relation

$|\psi(x_{1})+\psi(x_{2})|^{2}=|\psi(x_{1})|^{2}+|\psi(x_{2})|^{2}+\psi(x_{1})\psi^{*}(x_{2})+\psi(x_{2})\psi^{*}(x_{1}).$

The possibility of direct measurement of wavefunction is attributed to the
introduction of weak value by Aharonov, Albert and Vaidman (AAV) in 1988 AAV .
Lundeen et al first demonstrated experimentally the direct spatial
wavefunction reconstruction of photons via weak value in 2011 Lundeen , and
serious works on this topic has been reported since then wave1 ; wave2 ; wave3
; wave4 . In the original work of AAV, the weak value is defined based on the
framework of weak measurements. It has been realized recently that the weak
value is independent of measurement strength, which not only extent the
physical meaning of weak value, but also makes the practical measurement more
accurate wave5 . To date, however, direct measurement of superposition
wavefunction has not been reported.

Here we show that the direct measurement of superposition wavefunction based
on weak value is feasible, and the current experimental technologies is
already prepared for the target. With the ability to realize direct
measurement of superposition wavefunction, we get an important opportunity to
directly verify the Born rule.

##  II Weak value based direct superposition wavefunction measurement

The weak value of observable $\hat{A}$ is defined as

$<\hat{A}>_{w}=\frac{\langle\phi_{f}|\hat{A}|\psi_{i}\rangle}{\langle\phi_{f}|\psi_{i}\rangle},$
(1)

where $|\psi_{i}\rangle,|\phi_{f}\rangle$ represent the initial state and
post-selection state of the system respectively. If $\hat{A}$ is replaced with
position projector $|x\rangle\langle x|$ and $|\phi_{f}\rangle$ be zero
momentum state $|p=0\rangle$, then we obtain

$<|x\rangle\langle x|>_{w}=\frac{\langle p=0|x\rangle\langle
x|\psi_{i}\rangle}{\langle p=0|\psi_{i}\rangle}=C\cdot\psi_{i}(x)$ (2)

with constant $C=1/\langle p=0|\psi_{i}\rangle$. The real and imaginary part
of $<|x\rangle\langle x|>_{w}$ can be extracted directly by performing
suitable observables in experiment, which reveals the information of
wavefunction $\psi_{i}(x)$.

The extension to the case of superposition wavefunction is direct. Consider
$\hat{A}=|x_{1}\rangle\langle x_{1}|+|x_{2}\rangle\langle x_{2}|$, we have

$<|x_{1}\rangle\langle x_{1}|+|x_{2}\rangle\langle
x_{2}|>_{w}=C\cdot[\psi_{i}(x_{1})+\psi_{i}(x_{2})].$ (3)

The direct measurement of superposition wavefunction requires direct
measurement of superposition weak value $<|x_{1}\rangle\langle
x_{1}|+|x_{2}\rangle\langle x_{2}|>_{w}$. We now specify the physical
realization of the direct measurement of superposition weak value for more
general case. The interaction Hamiltonian between the system and the pointer
can be settled as

$\hat{H}=\theta\cdot\sum_{i}|x_{i}\rangle\langle
x_{i}|\otimes\hat{\sigma}_{y},$ (4)

where $\theta=\pi/2$ is the coupling strength and $\hat{\sigma}_{y}$ is the
Pauli operator. The unitary evolution corresponding to this Hamiltonian
$\hat{U}=e^{-i\hat{H}}$ implies that the qubit pointer be rotated $\pi$ along
$Y$ axis of the Bloch sphere if the system in the position $x_{i}$. Suppose
that the initial states of the system and pointer are $|\psi\rangle_{s}$ and
$|0\rangle_{p}$ respectively, the composite state after the interaction
becomes

$\begin{split}|\Psi\rangle_{sp}&=\hat{U}|\psi\rangle_{s}\otimes|0\rangle_{p}\\\
&=|\psi\rangle_{s}\otimes|0\rangle_{p}-\sqrt{2}\sum\psi_{s}(x_{i})|x_{i}\rangle\otimes|-\rangle_{p}\end{split}$
(5)

with $|-\rangle=(|0\rangle-|1\rangle)/\sqrt{2}$. The post-selection of the
system to the zero momentum state collapse the pointer into the state

$\begin{split}|\phi\rangle_{p}&=\langle p=0|\Psi\rangle_{sp}\\\
&=\phi_{s}(p=0)|0\rangle_{p}-\sqrt{2}\sum\psi_{s}(x_{i})\langle
p=0|x_{i}\rangle\otimes|-\rangle_{p}\\\
&=\frac{1}{C}[|0\rangle_{p}-\sqrt{2}<\sum|x_{i}\rangle\langle
x_{i}|>_{w}|-\rangle_{P}].\end{split}$ (6)

The real and imaginary part of superposition weak value
$<\sum|x_{i}\rangle\langle x_{i}|>_{w}$ can be directly extracted by the
proper measurements. Specifically, we have

$\begin{split}&\mathrm{Re}[<\sum|x_{i}\rangle\langle
x_{i}|>_{w}]=\frac{1}{2|C|^{2}}\langle\phi|(\hat{\sigma}_{x}+2|1\rangle\langle
1|)|\phi\rangle_{p}\\\ &\mathrm{Im}[<\sum|x_{i}\rangle\langle
x_{i}|>_{w}]=\frac{1}{2|C|^{2}}\langle\phi|\hat{\sigma}_{y}|\phi\rangle_{p}.\end{split}$
(7)

The above process for superposition wavefunction reconstruction is totally
same with the previous wavefunction reconstruction wave1 , but only with
$|x\rangle\langle x|$ be replaced with $\sum|x_{i}\rangle\langle x_{i}|$.

Figure 1: Experimental setup proposal for direct measurement pf superposition
wavefunction and thus directly verification of the Born rule. The single
photons are prepared in the polarization state
$|0\rangle_{p}=(|H\rangle+|V\rangle)/\sqrt{2}$ initially, and then reflected
by a phase-type SLM in which $\pi$ phase difference between $|H\rangle$ and
$|V\rangle$ can be added in the chosen areas. Photons that reflected from the
chosen areas are converted into the orthogonal polarization state
$|1\rangle_{p}=(|H\rangle-|V\rangle)/\sqrt{2}$, which realizes the position
$|x\rangle\langle x|$ measurements. The post-selection of zero momentum state
is completed by a Fourier lens and a pinhole. Finally the weak value is
directly extracted by measuring the proper observables using the polarization
analyser. The superposition wavefunction $\psi(x_{1})+\psi(x_{2})$ is obtained
by adding phases in areas $x_{1},x_{2}$ simultaneously, while
$\psi(x_{1}),\psi(x_{2})$ are obtained by adding phase only in $x_{1}$ and
$x_{2}$ respectively. PBS: polarization beam splitter, HWP: half wave plate,
QWP: quarter wave plate, SLM: spatial light modulator.

##  III Experimental proposal

The experimental setup to directly measure the superposition wavefunction and
thus verify the Born rule is almost the same as the one reported in our
previous work wave1 . There are, however, two differences. First, there is no
need to prepare specific spatial wavefunction of photons. Any state of the
system $|\psi\rangle_{s}$ is acceptable given that the preparation is the same
during all the runs. Second, there is no need to scan the position, and we
just add the $\pi$ phase in chosen area $x_{1}$ or $x_{2}$ or simultaneously
in sequence.

The setup, as shown in Fig. 1, consists of four parts, i.e., preparation of
initial states, measurement of position, post-selection of zero-momentum state
and extraction of weak value. Single photons can be generated by the
spontaneous parametric down-conversion process with the trigger Zhu . The
pointer is chosen as the polarization of photons, which is prepared in the
state $|0\rangle_{p}=(|H\rangle+|V\rangle)/\sqrt{2}$ by passing through a PBS
and a HWP rotated at $22.5^{\circ}$. The photons are reflected by a phase-only
SLM in which $\pi$ phase is added in the chosen areas. The polarization state
of photons that reflected from these areas will be converted into orthogonal
state $|1\rangle_{p}=(|H\rangle-|V\rangle)/\sqrt{2}$. This completes the
position $|x\rangle\langle x|$ measurement because we can know definitely
whether or not photons have passed the position $x$ according to the
polarization state of photons. The post-selection of zero momentum state is
realized by a Fourier lens and a pinhole placed in its focal plane. The weak
value is extracted by measuring the polarization pointer with a polarization
analyser consists of HWP, QWP and PBS. The direct quantites that we obtained
in experiment is the expectation values
$\langle\phi|\hat{\sigma}_{x}|\phi\rangle_{p},\langle\phi|\hat{\sigma}_{y}|\phi\rangle_{p}$
and $\langle\phi|1\rangle\langle 1|\phi\rangle_{p}$. It should be noted that
$\hat{\sigma}_{x},\hat{\sigma}_{y}$ are defined in the basis of
$\\{|0\rangle,|1\rangle\\}$. Since weak value is proportional to the
wavefunction and the constant $C$ is fixed in experiment for all runs, we
could omit the constant and write the superposition wavefunction with the
observed expectation values as

$\begin{split}&\mathrm{Re}[\sum\psi(x_{i})]=\langle\phi|(\hat{\sigma}_{x}+2|1\rangle\langle
1|)|\phi\rangle_{p}\\\
&\mathrm{Im}[\sum\psi(x_{i})]=\langle\phi|\hat{\sigma}_{y}|\phi\rangle_{p}.\end{split}$
(8)

##  IV Discussion and Conclusion

Currently, there seems no technical difficulty to complete the suggested
experiment. However, it would not be easy to obtain accurate results. The key
is to reduce the system error introduced by the imperfect instruments as small
as possible. The statistical error can be reduced significantly by the
cumulative counting. The accurate experimental results will tells us to what
extent the Born rule applies.

In conclusion, we have introduced the weak value based direct measurement of
superposition wavefunction and its potential to the direct verification of the
Born rule. The specific experimental setup of photonic system to fulfill the
target is proposed.

###### Acknowledgements.

Meng-Jun Hu acknowledges the fruitful discussions with Prof. Yong-Sheng Zhang.
Meng-Jun Hu especially thanks the supported from the Bejing Academy of Quantum
information Sciences.

## References

  * (1) Y. Aharonov, D. Z. Albert, and L. Vaidman, How the result of a measurement of a component of the spin of a spin-1/2 particle can turn out to be 100, Phys. Rev. Lett 60, 1351 (1988).
  * (2) J. S. Lundeen, B. Sutherland, A. Patel, C, Stewart, and C. Bamber, Direct measurement of the quantum wavefunction, Nature 474, 188-191 (2011).
  * (3) C. R. Zhang, M. J. Hu, Z. B. Hou, J. F. Tang, J. Zhu, G. Y. Xiang, C. F. Li, G. C. Guo, and Y. S. Zhang, Direct measurement of the two-dimensional spatial quantum wave function via strong measurements, Phys. Rev. A 101, 012119 (2020).
  * (4) Z. M. Shi, M. Mirhosseini, J. Margiewicz, M. Malik, F. Rivera, Z. Y. Zhu, and R. W. Boyd, Scan-free direct measurement of an extremely high-dimensional photonic state, Optica 2, 4 (2015).
  * (5) C. R. Zhang, M. J. Hu, G. Y. Xiang, Y. S. ZHang, C. F. Li, and G. C. Guo, Direct Strong Measurement of a High-Dimensional Quantum State, Chin. Phys. Lett. 37, 080301 (2020).
  * (6) Y. Zhou, J. Zhao, D. Hay, K. McGonagle, R. W. Boyd, and Z. Shi, Direct Tomography of High-Dimensional Density Matrices for General Quantum States of Photons, Phys. Rev. Lett. 127, 040402 (2021).
  * (7) G. Vallone and D. Dequal, Strong Measurements Give a Better Direct Measurement of the Quantum Wave Function, Phys. Rev. Lett. 116, 040502 (2016).
  * (8) J. Zhu, M. J. Hu, S. M. Cheng, M. J. W. Hall, C. F. Li, G. C. Guo, and Y. S. Zhang, Experimental verification of anisotropic invariance for three-qubit states, Phys. Rev. A 99, 040103(R) (2019).


# An efficient truncation scheme for Eulerian and total Lagrangian SPH methods

Zhentong Wang<EMAIL_ADDRESS>Chi Zhang<EMAIL_ADDRESS>Oskar J.
Haidn<EMAIL_ADDRESS>Xiangyu Hu<EMAIL_ADDRESS>TUM School of
Engineering and Design, Technical University of Munich, Garching, 85747,
Germany  Huawei Technologies Munich Research Center, Munich, Germany

###### Abstract

In smoothed particle hydrodynamics (SPH) method, the particle-based
approximations are implemented via kernel functions, and the evaluation of
performance involves two key criteria: numerical accuracy and computational
efficiency. In the SPH community, the Wendland kernel reigns as the prevailing
choice due to its commendable accuracy and reasonable computational
efficiency. Nevertheless, there exists an urgent need to enhance the
computational efficiency of numerical methods while upholding accuracy. In
this paper, we employ a truncation approach to limit the compact support of
the Wendland kernel to 1.6h. This decision is based on the observation that
particles within the range of 1.6h to 2h make negligible contributions,
practically approaching zero, to the SPH approximation. To address integration
errors stemming from the truncation, we incorporate the Laguerre-Gauss kernel
for particle relaxation due to the fact that this kernel has been demonstrated
to enable the attainment of particle distributions with reduced residue and
integration errors [1]. Furthermore, we introduce the kernel gradient
correction to rectify numerical errors from the SPH approximation of kernel
gradient and the truncated compact support. A comprehensive set of numerical
examples including fluid dynamics in Eulerian formulation and solid dynamics
in total Lagrangian formulation are tested and have demonstrated that
truncated and standard Wendland kernels enable achieve the same level accuracy
but the former significantly increase the computational efficiency.

###### keywords:

Truncated Wendland kernel , Computational efficiency , Accuracy , Particle
relaxation , Laguerre-Gauss kernel , Eulerian SPH , Total Lagrangian SPH

##  1 Introduction

Smoothed particle hydrodynamics (SPH) emerges as a meshless methodology,
initially conceptualized by Lucy [2], Gingold and Monaghan [3] and has found
extensive traction across domains, including fluid dynamics [3], solid
mechanics [4], and various scientific and industrial applications [5, 6, 7].
Given that SPH relies on particle-based approximation facilitated by Gaussian-
like kernel functions within compact support, the choice of the kernel
functions significantly affects both the accuracy and computational
efficiency. The accuracy attained within SPH simulations depends on numerical
errors which encompass the smoothing error denoted as $E_{s}$ and the
integration error represented by $E_{r}$ [8], where the smoothing error is
introduced in the process of replacing the Dirac delta function with a smooth
kernel function and the error order is defined by vanishing moments [8], and
the integration error is introduced in the process of discretizing the
continuous volume into the summation of the volumes of all neighboring
particles. Importantly, the cumulative truncation error can be expressed as
$E_{t}=E_{s}+E_{r}$. In terms of computational efficiency, the size of the
cut-off radius determines the number of neighbouring particles, and different
kernels with different cut-off sizes have a noticeable effect on the
computational effort.

Within the SPH community, kernels such as Wendland [9], Gauss [10], and
B-Spline [11] are widely recognized and acknowledged and generally possess the
characteristic of satisfying the first vanishing moment [8], ensuring the
second-order smoothing error $E_{s}=\mathcal{O}(h^{2})$. While Gauss kernel as
a nature choice for SPH having the advantage of good accuracy [12], the
compact support must be truncated to $3h$ during implementation to ensure a
sufficiently small integration error and the truncated compact support is
larger than the non-truncated compact support of other kernels, leading to
poor computational efficiency. Moreover, Gauss kernel with truncated compact
support, similar to B-spline kernel, leads to a randomly disordered particle
system, where the particles exhibit a tendency to cluster together (known as
the pair-instability problem). Thus, in recent years, Wendland kernel with the
compact support of $2h$ has become the mainstream due to its capacity to
overcome the pair-instability problem and maintain high computational
efficiency while achieving the same level of accuracy as the Gauss kernel. For
the given kernel, enhancing the computational efficiency in SPH can be
achieved by decreasing the resolution, i.e. reducing the total number of
particles in the whole computational domain, at the cost of decreasing the
accuracy. Therefore, the challenge lies in improving computational efficiency
without compromising the accuracy.

It is widely recognized that particle distribution has a significant impact on
the integration error [13, 1] and obtaining a high-quality particle
distribution becomes imperative and critical in preserving numerical accuracy.
In practical applications, particle relaxation has emerged as a remarkably
effective technique for generating body-fitted particles characterized by
complex geometries. While the utilization of Wendland kernel for particle
relaxation circumvents the pair-instability problem and resembles typical
liquid molecules in a microscopic description [8, 14], the quest for higher-
quality particle distribution leads to the introduction of Laguerre-Gauss
kernel investigated by Wang [1] as a replacement for Wendland kernel. In his
work, Wang [1] proves that the Laguerre-Gauss kernel has the property of
achieving a substantially smoother particle distribution, similar to the
lattice particle distribution, with significantly reduced relaxation residue
in the relaxation process. Besides, by analyzing the profile of Wendland
kernel function, it is evident that the value of the kernel function rapidly
diminishes as the distance from the center particle increases. At a distance
of $1.6h$, the value almost approaches zero, rendering the contribution of all
particles within the $1.6h$ to $2h$ range virtually negligible. Drawing
inspiration from this observation, we opt to exclusively utilize particles
within the truncated support domain to $1.6h$ for SPH approximation, leading
to the reduction of the total neighboring particles in this compact support
and a significant improvement in computational efficiency. Note that even
though the kernel values of particles within the support domains from $1.6h$
to $2h$ close to zero leading to little effect on the truncation error
intuitively, it is necessary to quantitatively analyze the impacts of
truncating this region on the truncation error.

In this paper, we conducts a quantitative analysis of truncation errors
employing both the standard Wendland kernel (with full compact support as
$2h$) and the truncated Wendland kernel (with compact support as $1.6h$). This
analysis is performed with and without the kernel gradient correction across
three distinct particle distributions, which encompass the lattice
distribution and the relaxed distributions utilizing both the standard
Wendland and the truncated Laguerre-Gauss kernels for particle relaxation.
Furthermore, the truncated Wendland kernel is implemented in a series of
numerical examples, spanning fluid dynamics within the Eulerian SPH framework
and solid dynamics in the total Lagrangian SPH framework to assess the
numerical accuracy and computational efficiency in comparison to the standard
Wendland kernel.

The structure of paper begins with an in-depth presentation of SPH
approximation and the characteristics of kernel in Chapter 2. In particular,
the truncated Wendland kernel is discussed in detail and subjected to a
quantitative analysis. Moving on to Chapter 3, the conservation equation and
its corresponding discretization in Eulerian SPH method are provided. In
addition, the Riemann solvers, dissipation limiters [15, 16, 17, 18] and the
kernel gradient correction [19] are introduced. The governing equations and
discretization forms in total Lagrangian SPH method are also explained. In
Chapter 4, a series of numerical examples are used to meticulously investigate
the numerical accuracy and computational efficiency of the truncated Wendland
kernel.

##  2 Truncated kernel for SPH method

###  2.1 SPH approximations

In SPH method, the gradient of a function $f(\mathbf{r})$ at the position
$\mathbf{r}_{i}$ can be approximated as

$\begin{split}\nabla f(\mathbf{r}_{i})\approx\int\limits_{\Omega}\nabla
f(\mathbf{r})W(\mathbf{r}_{i}-\mathbf{r},h)dV=-\int\limits_{\Omega}f(\mathbf{r})\nabla
W(\mathbf{r}_{i}-\mathbf{r},h)dV\approx-\sum_{j}f(\mathbf{r}_{j})\nabla
W_{ij}V_{j}.\end{split}$ (1)

Here, $V_{j}$ is the particle volume and
$W_{ij}=W(\mathbf{r}_{i}-\mathbf{r}_{j},h)$, and the gradient of kernel
$\nabla W_{ij}=\frac{\partial W_{ij}}{\partial r_{ij}}\mathbf{e}_{ij}$. Note
that the first approximation of Eq. (1) introduces the smoothing error $E_{s}$
due to the smooth kernel function replacing the Dirac delta function, and the
second approximation introduces the integration error $E_{r}$, and then the
truncation error $E_{t}$ with $E_{t}=E_{s}+E_{r}$. Following Ref.[20], a weak
form of $\nabla f(\mathbf{r}_{i})$ can be rewritten as

$\nabla f_{i}=\nabla f_{i}+f_{i}\nabla 1\approx-2\sum_{j}\bar{f}_{ij}\nabla
W_{ij}V_{j}$ (2)

to discretize the conservation equations due to its conservation property.
Here, $\bar{f}_{ij}=(f_{i}+f_{j})/2$ denotes the particle average value of
particles $i$ and $j$. Note that Eq. (2) is employed in subsequent sections
for the error estimation

###  2.2 Properties of kernel function

As discussed in Ref. [9, 21, 22, 23], a smooth kernel function must satisfy
several conditions which are three important properties:

$(1)$ Normalization condition, i.e unity condition, requires that the kernel
function satisfies

$\int\limits_{V}W(\mathbf{r},h)dV=1.$ (3)

$(2)$ When the smooth length tends to 0, the kernel function follows the Dirac
function property, that is

$\lim_{h\to 0}W(\mathbf{r},h)=\delta(\mathbf{r}).$ (4)

$(3)$ The compact condition requires

$W(\mathbf{r},h)=0,\text{when}\left|\mathbf{r}\right|\textgreater\kappa h,$
(5)

where $\kappa$ determine the size of compact support.

###  2.3 Truncated Wendland kernel

####  2.3.1 Truncated kernel basics

Following Ref.[9], Wendland kernel reads

$W(s,h)=\alpha_{d}\begin{cases}(1-s/2)^{4}(2s+1)&\text{ if }0\leq s\leq\kappa
h\\\ 0&\text{ if }s>\kappa h\end{cases},$ (6)

where $s=\left|\mathbf{r}\right|$ and $\alpha_{d}$ is the normalized
coefficient in $d$ dimensional space with the value of $\alpha_{1}=3/(4h)$,
$\alpha_{2}=7/(4\pi h^{2})$ and $\alpha_{3}=21/(16\pi h^{3})$. In the SPH
community, the constant parameter $\kappa$ is usually adopted as $2$, i.e. the
cut-off radius is $2h$. As shown in Figure 1, the kernel function curve (red
one) and its derivative (blue one) are plotted with the black and green
circles representing the compact support domain with cut-off size of $2h$ and
$1.6h$, respectively.

Figure 1: Wendland kernel: function (red) and its derivative (blue) in the
compact support and two different compact support size with $\kappa=2$ and
$\kappa=1.6$ are plotted as black and green circles, respectively.

It is straightforward to find that $W(s,h)$ is very close to $0$ in the range
of support domains where particle distance is larger than $1.5h$, implying
particles in this range contribute essentially nothing to the normalization
condition. Through this observation, we can truncate the compact support size
from $2h$ to $1.6h$, which leads to a reduction in the number of neighboring
particles, and therefore optimizes computational efficiency. Note that the
Wendland kernel with truncated compact support of $1.6h$ satisfies the
normalization condition due to the observation that little distribution of the
particles located in the region where particle distance is larger than $1.6h$,
which will be discussed quantitatively in section 2.3.3. While the first
derivative of the kernel is not very close to zero in the truncated region,
the kernel gradient correction matrix therefore introduced in section 2.3.2
eliminates and compensates the truncation error. It is worth noting that the
present trunction scheme does not affect another two conditions of Eqs. (4)
and (5).

####  2.3.2 Preparation of error analysis

In subsequent content, the error analysis is undertaken within the prescribed
particle distribution, encompassing both the lattice distribution and the
relaxed distributions. The particle relaxation technique [24] on the latter
distributions is defined by the following expression:

$\mathbf{a}_{p,i}=-\frac{2}{\rho}_{i}\sum_{j}p_{0}\nabla W_{ij}V_{j},$ (7)

with $p_{0}$ denoting the constant pressure. It is noteworthy that the
acceleration attains the value of zero when the particles attain their steady
states. Besides, the relaxed distributions in the error analysis contain the
particle relaxation using standard Wendland kernel and Laguerre-Gauss kernel
[1] which reads

$W^{LG}{(s,h)}=\alpha_{d}\begin{cases}(1-\frac{s^{2}}{2}+\frac{s^{4}}{6})e^{-s^{2}}&\text{
if }0\leq s\leq 2\\\ 0&\text{ if }s>2\end{cases},$ (8)

with $\alpha_{d}$ denoting the normalized coefficient where $d$ is dimensional
space and the value of $\alpha_{1}=8/(5\times\sqrt{\pi})$, $\alpha_{2}=3/\pi$
and $\alpha_{3}=8/\pi^{3/2}$. Moreover, to remedy the truncation error and the
approximation of kernel gradient, the kernel gradient correction matrix [19]

$\mathbf{B}_{i}=-\left(\sum_{j}\mathbf{r}_{ij}\otimes\nabla
W_{ij}V_{j}\right)^{-1}$ (9)

is therefore implemented to improve the numerical accuracy.

####  2.3.3 Error analysis

To analyze the truncation error, we introduce $L_{2}$ error [25] of the
particle approximation and its first derivative given by

$L_{2}(\phi)=\frac{(\sum\limits_{i}(\left|\phi_{i}^{\text{analytical}}-\phi_{i}^{\text{numerical}}\right|)^{2})^{\frac{1}{2}}}{(\sum\limits_{i}(\left|\phi_{i}^{\text{analytical}}\right|)^{2})^{\frac{1}{2}}},$
(10)

where $\phi$ represents the function $f$ or its derivative $df/dx$. Here,
$\phi^{\text{numerical}}$ is approximated as the weak form in Eq.(2). In the
case, we investigate a circular geometry characterized by a diameter of $D=2$,
with a smoothing length of $h=1.3dp$, where $dp$ represents the initial
particle spacing. Three different resolutions, namely $dp=D/10$, $dp=D/20$,
and $dp=D/40$, are employed for the convergence study. For the clarity
throughout our analysis, we consider three different particle distributions: a
lattice distribution, a relaxed distribution using the standard Wendland
kernel $(\kappa=2.0)$, and a relaxed distribution using the Laguerre-Gauss
kernel. We refer to the errors using the standard Wendland kernel, both with
and without the kernel gradient correction, as ”SW With B” and ”SW Without B”,
respectively. Similarly, the errors using the truncated Wendland kernel
$(\kappa=1.6)$ are referred to as ”TW With B” and ”TW Without B”, with and
without the kernel gradient correction.

We firstly explore the normalization condition by investigating the SPH
approximation for $\phi$ as a constant function $f(x)=1$ with the resolution
$dp=D/10$. Note that the summation of the kernel weight is independent on the
gradient operator, and therefore we only consider the ”SW without B” and ”TW
without B”. Table 1 lists the errors of SPH approximation for the constant
function $f(x)=1$ in different particle distributions with the resolution
$dp=D/10$. The results show the error of the truncated Wendland kernel is
almost same as that of the standard kernel, validating that the present
trunction scheme does not affect the normalization condition of Eqs. (3).

Table 1: Error analysis of SPH approximation for an constant function $\phi:f(x)=1$ in different particle distributions with the resolution $dp=D/10$. Particle distribution |  | $L_{2}(\phi)$
---|---|---
| SW without B |  | TW without B
Lattice |  | 0.010 |  | 0.006
Wendland-based |  | 0.034 |  | 0.028
Lagueree-Gauss-based |  | 0.037 |  | 0.031

As discussed in early studies in Refs. [9, 8] that Wendland kernel exhibits
2nd-order smoothing error, to investigate the smoothing error of the truncated
Wendland kernel, we expand function $f(\mathbf{r})$ using Taylor series
expansion and rewrite Eq. (1) as

$\displaystyle\int\limits_{\Omega}f(\mathbf{r})W(\mathbf{r}_{i}-\mathbf{r},h)dV=$
(11) $\displaystyle\quad
f\left(\mathbf{r}_{i}\right)\int\limits_{\Omega}W\left(\mathbf{r}-\mathbf{r}_{i},h\right)dV+\nabla
f\left(\mathbf{r}_{i}\right)\int\limits_{\Omega}\left(\mathbf{r}-\mathbf{r}_{i}\right)W\left(\mathbf{r}-\mathbf{r}_{i},h\right)dV+\mathcal{O}(h^{2}).$

With the normalization condition and even function condition [23] in hand, it
is easy to derive

$\displaystyle\int\limits_{\Omega}f(\mathbf{r})W(\mathbf{r}_{i}-\mathbf{r},h)dV$
$\displaystyle=f\left(\mathbf{r}_{i}\right)+\mathcal{O}(h^{2}),$ (12)

indicating that the truncated Wendland kernel remains 2nd-order smoothing
error.

Figure 2: The convergence study of error for kernel gradient on lattice
distribution (upper panel) and the relaxed distributions using Wendland
(middle panel) and Laguerre-Gauss (bottom panel) kernels, respectively.

Moreover, in the context of assessing convergence and comparing errors between
employing the standard Wendland kernel and the truncated Wendland kernel, we
test the SPH approximation for the derivative of the function $df/dx$ with the
exponential function given by

$f(x)=e^{-\frac{x^{2}}{0.1}},$ (13)

with the resolutions including $dp=D/10$,$D/20$ and $D/40$. Figure 2
illustrates the $L_{2}$ norm of the SPH approximation for the derivative of an
exponential function across three distinct particle distributions, including
the lattice distribution and relaxed distributions employing the Wendland and
Laguerre-Gauss kernels and the two lines with scopes equal to $1$ and $2$.

In the top panel of Figure 2, which pertains to the lattice distribution, it
becomes evident that the truncation scheme introduces the extra error by
comparing the lines labeled ”SW Without B” and ”TW Without B”. However, the
introduction of the kernel gradient correction effectively compensates for the
errors originating from SPH approximation and the truncated cut-off radius
observed in the lines of ”SW With B” and ”TW With B”. Consequently, the
results using the standard and truncated Wendland kernels with the kernel
gradient correction obtain the very similar error and approach nearly second-
order convergence.

However, in the middle panel of Figure 2, it is evident that, while the kernel
gradient correction still reduces errors, a clear trend emerges which is that,
as the resolution transitions from $dp=D/20$ to $dp=D/40$, the error using the
truncated Wendland kernel cannot be reduced further. It is noteworthy that,
within the given relaxed distribution using the standard Wendland kernel,
employing the truncation method results in significantly higher errors. This
is primarily attributed to the truncation scheme’s strong reliance on the
quality of the particle distribution, specifically, a lower relaxation
residue. Consequently, the relaxed distribution using the standard Wendland
kernel is unsuitable for the proposed truncation scheme.

In the bottom panel of Figure 2, the effect of the kernel gradient correction
on error reduction is evident within the given relaxed distribution using the
Laguerre-Gauss kernel. While the error when employing the truncated Wendland
kernel (as seen in the line ”TW Without B”) is larger than that when using the
standard Wendland kernel (as indicated in the line ”SW Without B”), the kernel
gradient correction introduced brings the errors for both the standard and
truncated Wendland kernels to the same level and achieves the second-order
convergence (as depicted in the lines ”SW With B” and ”TW With B”). This
result notably differs from the findings in the middle panel, signifying that
the particle distribution employing the Laguerre-Gauss kernel is best suited
for the truncation scheme.

##  3 Governing equations and SPH discretizations

###  3.1 Eulerian SPH for fluid dynamics

####  3.1.1 Governing equations

The Euler equation can be written as

$\frac{\partial\mathbf{U}}{\partial t}+\nabla\cdot\mathbf{F}(\mathbf{U})=0,$
(14)

where the conserved variables $\mathbf{U}$ and the corresponding fluxes
$\mathbf{F}(\mathbf{U})$ are described as

$\mathbf{U}=\left[\begin{array}[]{c}\rho\\\ \rho u\\\ \rho v\\\ \rho w\\\
E\end{array}\right]\quad\text{and}\quad\mathbf{F}=\left[\begin{array}[]{c}\rho
u\\\ \rho u^{2}+p\\\ \rho uv\\\ \rho uw\\\
u(E+p)\end{array}\right]+\left[\begin{array}[]{c}\rho v\\\ \rho vu\\\ \rho
v^{2}+p\\\ \rho vw\\\ v(E+p)\end{array}\right]+\left[\begin{array}[]{c}\rho
w\\\ \rho wu\\\ \rho wv\\\ \rho w^{2}+p\\\ w(E+p)\end{array}\right],$ (15)

where $u$, $v$ and $w$ are the components of velocity, $\rho$ and $p$ are the
density and pressure, respectively, and $E=\frac{\rho{\mathbf{v}}^{2}}{2}+\rho
e$ the total energy with $e$ denoting the internal energy. The equation is
closed with an equation of state (EOS). we apply the EOS of

$p=\rho(\gamma-1)e,$ (16)

with the heat capacity ratio $\gamma=1.4$ for compressible flows whose sound
speed is calculated by

$c^{2}=\frac{\gamma p}{\rho}.$ (17)

For incompressible flow, we apply the artificial EOS

$p=c^{2}(\rho-\rho_{0}).$ (18)

Here, following the weakly-compressible assumption, $\rho_{0}$ is the
reference density and the artificial sound speed $c=10U_{max}$ with $U_{max}$
denoting the maximum flow speed to limit the density variation less than 1%.
Note that the energy equation described in Eq. (15) is inactive during weakly-
compressible flow simulations.

####  3.1.2 Riemann-based Eulerian SPH discretization

Following Ref. [26], the Euler equation can be discretized as

$\left\\{\begin{array}[]{l}\frac{\partial}{\partial
t}\left(w_{i}\rho_{i}\right)+2w_{i}\sum_{j}w_{j}(\rho\mathbf{v})^{*}_{E,ij}\cdot\nabla
W_{ij}=0\\\ \frac{\partial}{\partial
t}\left(w_{i}\rho_{i}\mathbf{v}_{i}\right)+2w_{i}\sum_{j}w_{j}\left[(\rho\mathbf{v}\otimes\mathbf{v})^{*}_{E,ij}+p^{*}_{E,ij}\mathbb{I}\right]\cdot\nabla
W_{ij}=0\\\ \frac{\partial}{\partial
t}\left(w_{i}E_{i}\right)+2w_{i}\sum_{j}w_{j}\left[(E\mathbf{v})^{*}_{E,ij}+(p\mathbf{v})^{*}_{E,ij}\right]\cdot\nabla
W_{ij}=0\end{array},\right.$ (19)

where $w$, $\mathbf{v}$ and $\mathbb{I}$ are the volume of particle, velocity
and identity matrix, respectively, and terms $()^{*}_{E,ij}$ are solution of
the Riemann problem.

For the Riemann solution [15, 16, 17], the smallest and the largest wave
speeds are $S_{l}$ and $S_{r}$, respectively, and the middle speed $S_{\ast}$
is the discontinuity separating the solution into
$(\rho_{l}^{\ast},u_{l}^{\ast},p_{l}^{\ast})$ and
$(\rho_{r}^{\ast},u_{r}^{\ast},p_{r}^{\ast})$. For compressible flows, we
apply the HLLC Riemann solver with $S_{l}$ and $S_{r}$ given by

$\left\\{\begin{array}[]{l}S_{l}=u_{l}-c_{l}\\\
S_{r}=u_{r}+c_{r}\end{array},\right.$ (20)

the middle wave speed

$S_{\ast}=\frac{\rho_{r}u_{r}\left(S_{r}-u_{r}\right)+\rho_{l}u_{l}\left(u_{l}-S_{l}\right)+p_{l}-p_{r}}{\rho_{r}\left(S_{r}-u_{r}\right)+\rho_{l}\left(u_{l}-S_{l}\right)}.$
(21)

Then, the conserved variables in the intermediate region can be obtained as

$\left\\{\begin{array}[]{l}p^{*}=p_{l}+\rho_{l}\left(u_{l}-S_{l}\right)\left(u_{l}-u^{*}\right)=p_{r}+\rho_{r}\left(S_{r}-u_{r}\right)\left(u^{*}-u_{r}\right)\\\
\mathbf{v}_{l/r}^{*}=u^{*}\mathbf{e}_{ij}+\left[\frac{1}{2}(\mathbf{v}_{l}+\mathbf{v}_{r})-\frac{1}{2}({u}_{l}+{u}_{r})\mathbf{e}_{ij}\right]\\\
\rho_{l/r}^{*}=\rho_{l/r}\frac{\left(S_{l/r}-\left|\mathbf{v}_{l/r}\right|\right)}{\left(S_{l/r}-u^{*}\right)}\\\
E_{l/r}^{*}=\frac{\left(S_{l/r}-\left|\mathbf{v}_{l/r}\right|\right)E_{l/r}-p_{l/r}\left|\mathbf{v}_{l/r}\right|+p^{*}u^{*}}{S_{l/r}-u^{*}}\end{array},\right.$
(22)

with ${u}^{*}={S}_{\ast}$. For incompressible flows, we apply the linearised
Riemann solver [17] with

$\left\\{\begin{array}[]{l}u^{*}=\frac{u_{l}+u_{r}}{2}+\frac{1}{2}\frac{\left(p_{l}-p_{r}\right)}{\bar{\rho}\bar{c}}\\\
p^{*}=\frac{p_{l}+p_{r}}{2}+\frac{1}{2}\bar{\rho}\bar{c}\left(u_{l}-u_{r}\right)\end{array},\right.$
(23)

where $\bar{\rho}=\frac{1}{2}(\rho_{l}+\rho_{r})$ and
$\bar{c}=\frac{1}{2}(c_{l}+c_{r})$

To decrease the numerical dissipation, following Refs. [7, 18], the HLLC
Riemann solver of Eq. (22) adopts the dissipation limiter as

$\left\\{\begin{array}[]{l}u^{*}=\frac{\rho_{l}u_{l}c_{l}+\rho_{r}u_{r}c_{r}}{\rho_{l}c_{l}+\rho_{r}c_{r}}+\frac{p_{l}-p_{r}}{\rho_{l}c_{l}+\rho_{r}c_{r}}\beta^{2}_{HLLC}\\\
p^{*}=\frac{p_{l}+p_{r}}{2}+\frac{1}{2}\beta_{HLLC}\left[\rho_{r}c_{r}\left(u^{*}-u_{r}\right)-\rho_{l}c_{l}\left(u_{l}-u^{*}\right)\right]\end{array}.\right.$
(24)

Here, the limiter $\beta_{HLLC}$ reads

$\beta_{HLLC}=\min\left(\upeta_{HLLC}\max(\frac{u_{l}-u_{r}}{\bar{c}},0),1\right),$
(25)

with $\upeta_{HLLC}=1$ according to experimental and numerical examples. For
the linearised Riemann solver of Eq. (23), the dissipation limiter is
introduced as

$\left\\{\begin{array}[]{l}u^{*}=\frac{u_{l}+u_{r}}{2}+\frac{1}{2}\beta^{2}_{linearisd}\frac{\left(p_{L}-p_{R}\right)}{\bar{\rho}\bar{c}}\\\
p^{*}=\frac{p_{l}+p_{r}}{2}+\frac{1}{2}\beta_{linearisd}\bar{\rho}\bar{c}\left(u_{l}-u_{r}\right)\end{array}.\right.$
(26)

where the limiter $\beta_{linearisd}$ is defined as

$\beta_{linearisd}=\min\left(\upeta_{linearisd}\max(\frac{u_{l}-u_{r}}{\bar{c}},0),1\right),$
(27)

with $\upeta_{linearisd}=15$.

To compensate the errors introduced from SPH approximation for the kernel
gradient and truncated compact support in this paper, the kernel gradient
correction in Eq. (9) is introduced to rewrite the kernel gradient as

${\nabla}^{{}^{\prime}}W_{ij}=\frac{\mathbf{B}_{i}+\mathbf{B}_{j}}{2}\nabla
W_{ij}$ (28)

to replace original $\nabla W_{ij}$ of Eq. (19).

###  3.2 Total Lagrangian SPH for solid dynamics

For solid dynamics, the conservation equations in total Lagrangian framework
read

$\left\\{\begin{array}[]{l}\rho=J^{-1}\rho^{0}\\\
\rho^{0}\frac{d\boldsymbol{v}_{i}}{dt}=\nabla^{0}\cdot\mathbb{P}^{\mathrm{T}}\end{array},\right.$
(29)

with $\rho$, $\mathbb{P}$ and $T$ are the density, the first Piola-Kirchhoff
stress tensor and the matrix transposition operator, respectively. The
superscript $0$ denote the value in underformed reference. Also,
$J=det(\mathbb{F})$ with $\mathbb{F}$ being the deformation gradient tensor.
Then, The first Piola-Kirchhoff stress tensor
$\mathbb{P}=\mathbb{F}\mathbb{S}$ with $\mathbb{S}$ denoting the second Piola-
Kirchhoff stress tensor given by [27]

$\displaystyle\mathbb{S}=\lambda\operatorname{tr}(\mathbb{E})\mathbb{I}+2G\mathbb{E}.$
(30)

Here, $\lambda$ and $G$ the Lam$\grave{e}$ parameter and the shear modulus,
respectively.

Following Ref. [27], the momentum equation of Eq. (29) can be discretized in
total Lagrangian SPH form as

$\rho_{i}^{0}\frac{d\boldsymbol{v}_{i}}{dt}=\sum_{j}\left(\mathbb{P}_{i}\mathbf{B}_{i}^{0^{\mathrm{T}}}+\mathbb{P}_{j}\mathbf{B}_{j}^{0^{\mathrm{T}}}\right)\cdot\nabla^{0}W_{ij}V_{j}^{0},$
(31)

where $\nabla^{0}W_{ij}=\frac{\partial
W\left(\boldsymbol{r}_{ij}^{0},h\right)}{\partial{r}_{ij}^{0}}\boldsymbol{e}_{ij}^{0}$
and $\mathbf{B}^{0}$ present the kernel gradient and kernel correction matrix
in the initial configuration, respectively. It is worth noting that the
deformation tensor can be updated by using its change rate defined as

$\frac{d\mathbb{F}}{dt}=\sum_{j}V_{j}^{0}\left(\boldsymbol{v}_{j}-\boldsymbol{v}_{i}\right)\cdot\nabla^{0}W_{ij}\mathbf{B}_{i}^{0}.$
(32)

##  4 Numerical results

In this section, we test a set of numerical cases including fluid and solid
dynamics to verify the numerical accuracy and computational efficiency of the
proposed truncated Wendland kernel. Note that the smoothing length is applied
as $h=1.3dp$ in fluid dynamics and $1.15dp$ in solid dynamics. For clarity,
the Eulerian and total Lagrangian methods using standard Wendland kernel with
cut-off radius as $2h$ (denoted as SW) and truncated Wendland kernel with cut-
off radius as $1.6h$ (denoted as TW), respectively. In all cases, the
computations are all performed on an Intel Core i7-10700 2.90 GHz 8-core
desktop computer and the CPU wall-clock times required by standard and
truncated Wendland kernel denoted as $\text{T}_{\text{SW}}$ and
$\text{T}_{\text{TW}}$, respectively. Here, we establish the parameter of
computational time rate, denoted as $\alpha$, which is given by the ratio of
the CPU wall-clock times using TW to those using SW, i.e. $\alpha={(\text{CPU
wall-clock times})_{\text{TW}}}/{(\text{CPU wall-clock times})_{\text{SW}}}$.
A lower value of this parameter indicates enhanced computational efficiency
when utilizing the TW.

###  4.1 Double Mach reflection of a strong shock

In this section, we undertake a two-dimensional investigation, specifically
the double Mach reflection phenomenon resulting from a high-strength
shockwave. This examination serves the purpose of a correctness analysis by
using SW and TW in the context of compressible flow dynamics.

Figure 3: Double Mach reflection of a strong shock: Density contour and its
zoom-in view ranging from $1.3$ to $24.0$ obtained by SW (top panel) and TW
(bottom panel) with the resolution $dp=1/125$.

Following Ref. [28], we define our computational domain as
$(x,y)\in[0,4]\times[0,1]$ and the initial conditions are prescribed as
follows

$(\rho,u,v,p)=\begin{cases}(1.4,0,0,1)&y\leq 1.732(x-0.1667)\\\
(8,7.145,-4.125,116.8333)&\text{ otherwise }\end{cases}.$ (33)

The simulation is performed until the final time $t=0.2$. Furthermore, there
exists an initial configuration wherein a Mach 10 shockwave propagates to the
right, its starting point being $(x,y)=(0.1667,0)$, and it maintains a
60-degree inclination relative to the x-axis. The computational domain is
bounded by distinct conditions: $(1)$ The bottom boundary, spanning from
$x=1/6$ to $x=4$, is defined as a reflective wall boundary; $(2)$ The left-
hand boundary is governed by the post-shock boundary condition; $(3)$ The
right boundary at $x=4$ adheres to a zero-gradient condition. Besides, the
spatial resolution is applied as $dp=1/125$.

Figure 3 shows the density contour ranging from $1.3$ to $24.0$ obtained by SW
and TW with the resolution $dp=1/125$ at the finial time $t=0.2$, indicating
that both methods successfully capture key flow features, such as the Mach
stem and the near-wall jet and obtain equivalent accuracy of the results.
Moreover, the CPU wall-clock times demanded by SW and TW throughout the entire
process with the resolution $dp=1/125$ are $226.20s$ and $137.27s$,
respectively. This indicates that the computational time using the latter is
approximately 60% that of the former, highlighting the significantly improved
computational efficiency of the latter over the former.

###  4.2 Double shear layer flow

In this section, the simulation involves a double periodic shear layer in an
inviscid flow $(\mu=0)$ across a square domain with unit edge length. The
initial flow depicts a shear layer with finite thickness and a minor vertical
velocity perturbation and the initial velocity field [29, 7] is defined as
follows:

$u_{x}=\begin{cases}\tanh((y-0.25)/\rho),&\text{for }y\leq 0.5,\\\
\tanh((0.75-y)/\rho),&\text{for }y>0.5,\end{cases}\quad\text{and}\quad
u_{y}=\delta\sin(2\pi x),$ (34)

where $\rho=1/30$ and $\delta=0.05$.

Figure 4: Double shear layer flow : The vorticity contour ranging from $-25$
to $25$ obtained by SW (left panel) and TW (right panel) with the spatial
resolution $dp=1/128$ at different times $t=0.0$, $0.8$ and $1.2$.

Figure 4 presents the vorticity contour ranging from $-25$ to $25$ obtained by
SW and TW with $dp=1/128$ at different times $t=0.0$, $0.8$ and $1.2$. It can
be noted that the truncated kernel is able to obtain results as smooth as the
standard kernel, and the results are highly similar, which proves that kernel
with or without truncation can obtain the same accuracy. Besides, the CPU
wall-clock times required by SW and TW in the whole process under the
resolution $dp=1/128$ are $29.99s$ and $17.83s$, respectively, and then
$\alpha$ is calculated as $0.59$, illustrating that the considerably improved
computational efficiency of using the proposed kernel.

###  4.3 Lid-driven cavity flows with two shapes

In this section, we delve into two-dimensional lid-driven cavity flows with
varying geometries to conduct a computational efficiency comparison between SW
and TW. Initially, we employ a straightforward square cavity, with the
computed results subjected to scrutiny against the benchmark results presented
by Ghia et al. [30] to validate their accuracy. Furthermore, we assess the
capability of the proposed kernel in handling intricate geometries by
conducting simulations on a semi-circular cavity. The outcomes of this
simulation are compared with the reference results obtained by Glowinski et
al. [31]. Here, in both situations, the simulations are run until the final
time $t=30$ and the upper moving wall is assigned a constant velocity
$U_{wall}=1.0$, while all other boundaries are maintained as non-slip wall
conditions.

####  4.3.1 Lid-driven square cavity flow

For the square cavity, the geometry and boundary condtions are shown in Figure
5 and three different spatial resolutions $dp=1/33$, $1/65$ and $1/129$ are
adopted for the convergence study. Also, the Reynolds number $Re=400$ is
applied in this case. From the Figure 5 (right panel), it can be noted that
the proposed kernel enables obtain the smooth results which are in good
agreement with the reference [30].

Figure 5: Lid-driven square cavity flow: The square geometry and boundary conditions. Besides, the velocity contour ranging from $2.2\times 10^{-5}$ to $0.96$ obtained by TW with the resolution as $dp=1/1/129$ under the Reynolds number $Re=400$. Figure 6: Lid-driven square cavity flow ($Re=400$): The horizontal velocity component $u$ along $x=0.5L$ (left panel) and the vertical velocity component $v$ along $y=0.5L$ (right panel) using SW with the resolution of $dp=1/33$ and TW with the resolutions of $dp=1/33$, $1/65$ and $1/129$, and the comparisons with the reference by Ghia [30]. Table 2: Lid-driven square cavity flow ($Re=400$): Comparison of computational efficiency using SW and TW with different spatial resolutions. Resolution | Kernel type | CPU wall-clock time (s) | $\alpha$
---|---|---|---
1/33 | SW | 6.36 | 0.76
TW | 4.81
1/65 | SW | 44.12 | 0.68
TW | 30.00
1/129 | SW | 371.78 | 0.60
TW | 222.71

Besides, Figure 6 presents the horizontal velocity component $u$ along
$x=0.5L$ and the vertical velocity component $v$ along $y=0.5L$ using SW and
TW with the different resolutions, and the comparisons with the reference by
Ghia [30]. It can be observed that employing TW allows to achieve the same
level of accuracy as SW under the resolution of $dp=1/33$, indicating that the
truncation method has a negligible impact on numerical accuracy. Also, the
results of using the proposed kernel achieve the second-order convergence with
the increase of the spatial resolutions. Furthermore, the CPU wall-clock times
using SW and TW required in the whole process are listed in Table 2, clearly
verifying that the truncation method has significant improvement in the
computational efficiency.

####  4.3.2 Lid-driven semi-circular cavity problem

In this section, we examine more intricate semi-circular geometries with a
diameter of $1$ to validate the efficacy of the proposed kernel in obtaining
precise and smooth results in flows characterized by complex geometries and
the resolution is applied as $dp=1/129$ to discretize the computational domain
as well as the Reynolds number $Re=1000$ is adopted.

Figure 7: Lid-driven semi-circular cavity flow: The semi-circular geometry and
boundary conditions. Besides, the velocity contour ranging from $7.4\times
10^{-5}$ to $0.96$ obtained by TW with the resolution as $dp=1/1/129$ under
the Reynolds number $Re=1000$. Figure 8: Lid-driven semi-circular cavity flow
($Re=1000$): The horizontal velocity component $u$ along $x=0.5L$ (left panel)
and the vertical velocity component $v$ along $y=-0.25L$ (right panel) using
SW and TW with the resolution of $dp=1/129$, and the comparisons with the
reference by Glowinski et al. [31].

Figure 7 presents the semi-circular geometry with the radius as $0.5$ and
boundary conditions and the velocity contour ranging from $7.4\times 10^{-5}$
to $0.96$ using TW with the resolution as $dp=1/1/129$ under the Reynolds
number $Re=1000$, showing that the velocity contour obtained by TW is smooth.
Also, Figure 8 illustrates the horizontal velocity component $u$ along
$x=0.5L$ and the vertical velocity component $v$ along $y=-0.25L$ using SW and
TW with the resolution as $dp=1/129$, and the comparisons with the reference
by Glowinski et al. [31], validating that truncated or standard kernels can
obtain the same level of the accuracy in the high resolution. Moreover, the
CPU wall-clock times simulated by SW and TW for the whole process with the
resolution $dp=1/128$ are $144.10s$ and $89.14s$, respectively, and further
$\alpha=0.64$, proving that the using the proposed kernel can significantly
improve computational efficiency.

###  4.4 Flow past a circular cylinder

In the section, the flow around a circular cylinder is considered to test
fluid-solid interaction with complex geometry. To verify the simulation
quantitatively, the drag and lift coefficient are given by

$C_{D}=\frac{2F_{D}}{\rho_{\infty}u_{\infty}^{2}A}\quad\text{and}\quad
C_{L}=\frac{2F_{L}}{\rho_{\infty}u_{\infty}^{2}A},$ (35)

where $F_{D}$ and $F_{L}$ are the drag and lift force respectively and $A$ is
the size of the cross-sectional area of the cylinder. In the case, the
computational domain is $[25D,15D]$ with the center of the cylinder located in
$(7.5D,7.5D)$ where the cylinder diameter $D=2$, Strouhal number
$St=fD/u_{\infty}$ and the Reynolds number $Re=\rho_{\infty}u_{\infty}D/\mu$
set as $100$. Besides, all boundaries are set as far-field boundary conditions
and we apply three different resolutions $dp=1/10$, $1/20$ and $1/30$ to study
its convergence.

Figure 9: Flow around a cylinder: Velocity contour ranging from $1.8\times 10^{-3}$ to $1.5$ using SW (left panel) and TW (right panel) kernels with $Re=100$ under the resolution $dp=1/10$ at $t=300$. Figure 10: Flow around a cylinder: Drag (left panel) and lift (right panel) coefficients using truncated Wendland kernel with $Re=100$. Table 3: Flow around a cylinder: Drag and lift coefficients from different experimental and numerical results at Re=100.  | Parameters
---
| $C_{D}$
---
| $C_{L}$
---
| $S_{t}$
---
White[32] | 1.46 | - | -
Khademinejad et al.[33] | 1.30 $\pm$ - |  $\pm$0.292 | 0.152
Chiu et al.[34] | 1.35 $\pm$ 0.012 |  $\pm$0.303 | 0.166
Le et al.[35] | 1.37 $\pm$ 0.009 |  $\pm$0.323 | 0.160
Brehm et al.[36] | 1.32 $\pm$ 0.010 |  $\pm$0.320 | 0.165
Liu et al.[37] | 1.35 $\pm$ 0.012 |  $\pm$0.339 | 0.165
TW | 1.34 $\pm$ 0.010 |  $\pm$0.340 | 0.172
Table 4: Flow around a cylinder: Comparison of computational efficiency using SW and TW with different spatial resolutions. Resolution | Kernel type | CPU wall-clock time (s) | $\alpha$
---|---|---|---
1/10 | SW | 2672.98 | 0.73
TW | 1949.30
1/20 | SW | 21397.50 | 0.72
TW | 15412.23
1/30 | SW | 68340.33 | 0.72
TW | 49256.36

Figure 10 presents the drag and lift coefficients using truncated Wendland
kernel with the resolutions $dp=1/10$, $1/20$ and $1/30$ under the Reynolds
number $Re=100$, showing that the drag coefficients attain a stable value
after an initial period of fluctuation, while the lift coefficient continues
to oscillate around zero. The converged result is documented in Table 3,
alongside various other experimental and numerical findings [32, 33, 34, 35,
36, 37] and the comparisons demonstrates the well agreement between the
converged result obtained by proposed method and the referenced data,
indicating its reasonability. Furthermore, the CPU wall-clock times required
by SW and TW in whole process with different resolutions are listed in Table
4, clearly showing that he computational time using truncated Wendland kernel
is approximately 70% of that using the standard Wendland kernel and the
computational efficiency improves significantly by using the former.

###  4.5 Twisting column

In this section, the 3D simulation twisting column is further studied to
verify the correctness of the proposed kernel. Following Refs. [38, 39, 40],
the column with the height $H=1m$ and length $L=6m$ is assumed to consist of a
nearly incompressible material, characterized by a density of
$\rho_{0}=1100Kg/m^{3}$, Young’s modulus and Poisson’s ratio are $E=17MPa$ and
$\nu=0.499$, respectively. Besides, the initial condition for the twisting
motion involves a sinusoidal rotational velocity field represented as
$\mathbf{\omega}=[0,\Omega_{0}\sin{(\pi y_{0}/2L)},0]$ with
$\omega_{0}=105rad/s$.

Figure 11: Twisting column: Mises stress ranging from $0MPa$ to $8.3MPa$ at
different time instants obtained by SW (top panel) and TW (bottom panel) with
initial rotational velocity $\mathbf{\omega}=[0,\Omega_{0}\sin{(\pi
y_{0}/2L)},0]$ with $\omega_{0}=105$ rad/s under the spatial resolution
$H/dp=10$ where $dp$ is the initial particle spacing. Figure 12: Twisting
column: Mises stress ranging from $0MPa$ to $9.8MPa$ at different time
instants obtained by TW with initial rotational velocity
$\mathbf{\omega}=[0,\Omega_{0}\sin{(\pi y_{0}/2L)},0]$ with $\omega_{0}=105$
rad/s under the different spatial resolutions as $H/dp=10$ (left panel),
$H/dp=20$ (middle panel) and $H/dp=30$ (right panel) where $dp$ is the initial
particle spacing at the time $t=0.0909s$.

Figure 11 portrays the Mises stress distribution of twisting column at
different timings using SW and TW with $\omega_{0}=105$ rad/s under the
resolution $H/dp=10$, showing that both kernels can obtain the smooth stress
distribution which are agreement well with the reference [40]. Furthermore,
three different resolutions are applied to verify its numerical stability
shown in Figure 12 and it can be observed that the smooth results can be
obtained by TW with the increase of the resolutions and seem as stable. To
compare the computational efficiency between thses both kernels, the CPU wall-
clock times required in the whole process are listed in Table 5, indicating
that the considerably improved computational efficiency can be achieved by the
proposed kernel.

Table 5: Twisting column: Comparison of computational efficiency using SW and TW with different spatial resolutions. Resolution | Kernel type | CPU wall-clock time (s) | $\alpha$
---|---|---|---
1/10 | SW | 27.55 | 0.67
TW | 18.52
1/20 | SW | 366.91 | 0.68
TW | 247.94
1/30 | SW | 1825.15 | 0.64
TW | 1173.56

###  4.6 Bending column

In this section, a $3D$ bending column is considered to test the accuracy and
computational efficiency of the proposed method. Following Refs. [41, 40], in
this case, we employ a rubber-like material with a length of $L=6m$ and a
height of $H=1m$ and the material’s cross-sectional area is square in shape
and remains fixed at the bottom during the simulation. Besides, the initial
velocity is $\mathbf{v}_{0}=10(\frac{\sqrt{3}}{2},\frac{1}{2},0)^{T}$, density
$\rho_{0}=1100kg/m^{3}$, Young’s modulus $E=17MPa$ and Poisson’s ratio
$\nu=0.45$. For the convergence study, we apply three different resolutions
$dp=1/6$, $1/12$ and $1/24$.

Figure 13: Bending column: The Mises stress from $0MPa$ to $2.8MPa$ using TW with the spatial resolution $H/dp=12$ at the different instants. Figure 14: Bending column: The vertical position z located at the top of the column [40] using SW and TW with three different resolutions and the comparisons with the reference data from Aguirre et al. [42]. Table 6: Bending column: Comparison of computational efficiency using SW and TW with different spatial resolutions. Resolution | Kernel type | CPU wall-clock time (s) | $\alpha$
---|---|---|---
1/6 | SW | 1.45 | 0.69
TW | 1.01
1/12 | SW | 27.64 | 0.47
TW | 13.12
1/24 | SW | 626.31 | 0.56
TW | 352.11

Figure 13 presents the Mises stress from $0MPa$ to $2.8MPa$ using TW with the
spatial resolution $H/dp=12$ at the different timings, clearly showing that
the smooth results can be obtained by the proposed kernel well. Besides,
Figure 14 depicts the vertical position $z$ located at the top of the column
[40] using SW and TW with three different resolutions and the comparisons with
the reference data from Aguirre et al. [42]. From the results, the results
using two methods at the resolution $dp=1/6$ are almost coincide, verifying
that these methods enable obtain the same accuracy. Also, with the increase of
the spatial resolutions, the results converge rapidly and achieve second order
convergence. More importantly, the computational efficiency using SW and TW
are listed in Table 6, showing that the computational time using the latter is
almost half of that using the former, i.e. significantly improved
computational efficiency using truncated Wendland kernel.

##  5 Summary and conclusion

In this paper, the truncation error of the Wendland kernel, with its compact
support truncated to $1.6h$, has been subject to comprehensive quantitative
analysis, and the results have demonstrated that the error when using the
truncated Wendland kernel is comparable to that of the standard Wendland
kernel in conjunction with the kernel gradient correction. A series of
numerical tests encompassing both fluid and solid dynamics have been
conducted, and clearly show that the proposed truncation scheme results in
significantly improved computational efficiency without compromising numerical
accuracy. Note that the extent of this improvement does vary depending on the
specific scenario. Overall, the computational time required for calculations
involving the truncated Wendland kernel has been reduced by a minimum of 30%
compared to the demands of the standard Wendland kernel.

## References

  * [1] Z. Wang, B. Zhang, O. J. Haidn, X. Hu, A fourth-order kernel for improving numerical accuracy and stability in eulerian and total lagrangian sph, arXiv preprint arXiv:2309.01581 (2023).
  * [2] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The astronomical journal 82 (1977) 1013–1024.
  * [3] J. J. Monaghan, Simulating free surface flows with sph, Journal of computational physics 110 (2) (1994) 399–406.
  * [4] L. D. Libersky, A. G. Petschek, T. C. Carney, J. R. Hipp, F. A. Allahdadi, High strain lagrangian hydrodynamics: a three-dimensional sph code for dynamic material response, Journal of computational physics 109 (1) (1993) 67–75.
  * [5] P. Randles, L. D. Libersky, Smoothed particle hydrodynamics: some recent improvements and applications, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 375–408.
  * [6] S. M. Longshaw, B. D. Rogers, Automotive fuel cell sloshing under temporally and spatially varying high acceleration using gpu-based smoothed particle hydrodynamics (sph), Advances in Engineering Software 83 (2015) 31–44.
  * [7] Z. Wang, C. Zhang, O. J. Haidn, X. Hu, An eulerian sph method with weno reconstruction for compressible and incompressible flows, Journal of Hydrodynamics (2023) 1–12.
  * [8] S. Litvinov, X. Hu, N. A. Adams, Towards consistence and convergence of conservative sph approximations, Journal of Computational Physics 301 (2015) 394–401.
  * [9] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational Mathematics 4 (1995) 389–396.
  * [10] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Monthly notices of the royal astronomical society 181 (3) (1977) 375–389.
  * [11] G. R. Johnson, R. A. Stryk, S. R. Beissel, Sph for high velocity impact computations, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 347–373.
  * [12] J. P. Morris, Analysis of smoothed particle hydrodynamics with applications, Monash University Australia, 1996.
  * [13] N. J. Quinlan, M. Basa, M. Lastiwka, Truncation error in mesh-free particle methods, International Journal for Numerical Methods in Engineering 66 (13) (2006) 2064–2085.
  * [14] J. G. Kirkwood, E. M. Boggs, The radial distribution function in liquids, The Journal of Chemical Physics 10 (6) (1942) 394–402.
  * [15] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the hll-riemann solver, Shock waves 4 (1) (1994) 25–34.
  * [16] E. F. Toro, The hllc riemann solver, Shock waves 29 (8) (2019) 1065–1082.
  * [17] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science & Business Media, 2013.
  * [18] Z. Wang, C. Zhang, O. J. Haidn, N. A. Adams, X. Hu, Extended eulerian sph and its realization of fvm, arXiv preprint arXiv:2309.01596 (2023).
  * [19] V. Zago, L. J. Schulze, G. Bilotta, N. Almashan, R. Dalrymple, Overcoming excessive numerical dissipation in sph modeling of water waves, Coastal Engineering 170 (2021) 104018.
  * [20] C. Zhang, Y.-j. Zhu, D. Wu, N. A. Adams, X. Hu, Smoothed particle hydrodynamics: Methodology development and recent achievement, Journal of Hydrodynamics 34 (5) (2022) 767–805.
  * [21] M. Liu, G. Liu, Smoothed particle hydrodynamics (sph): an overview and recent developments, Archives of computational methods in engineering 17 (2010) 25–76.
  * [22] J. J. Monaghan, Smoothed particle hydrodynamics, Annual review of astronomy and astrophysics 30 (1) (1992) 543–574.
  * [23] X. Yang, S. Peng, M. Liu, A new kernel function for sph with applications to free surface flows, Applied Mathematical Modelling 38 (15-16) (2014) 3822–3833.
  * [24] Y. Zhu, C. Zhang, Y. Yu, X. Hu, A cad-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction, Journal of Hydrodynamics 33 (2) (2021) 195–206.
  * [25] A. Vergnaud, G. Oger, D. Le Touzé, Investigations on a high order sph scheme using weno reconstruction, Journal of Computational Physics 477 (2023) 111889\.
  * [26] J. Vila, On particle weighted methods and smooth particle hydrodynamics, Mathematical models and methods in applied sciences 9 (02) (1999) 161–209.
  * [27] C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, J. Wang, X. Hu, Sphinxsys: An open-source multi-physics and multi-resolution library based on smoothed particle hydrodynamics, Computer Physics Communications 267 (2021) 108066\.
  * [28] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, Journal of computational physics 54 (1) (1984) 115–173.
  * [29] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible navier-stokes equations, Journal of computational physics 85 (2) (1989) 257–283.
  * [30] U. Ghia, K. N. Ghia, C. Shin, High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method, Journal of computational physics 48 (3) (1982) 387–411.
  * [31] R. Glowinski, G. Guidoboni, T.-W. Pan, Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity, Journal of Computational Physics 216 (1) (2006) 76–91.
  * [32] F. M. White, J. Majdalani, Viscous fluid flow, Vol. 3, McGraw-Hill New York, 2006\.
  * [33] T. Khademinezhad, P. Talebizadeh, H. Rahimzadeh, Numerical study of unsteady flow around a square cylinder in compare with circular cylinder, in: Conference Paper, 2015.
  * [34] P.-H. Chiu, R.-K. Lin, T. W. Sheu, A differentially interpolated direct forcing immersed boundary method for predicting incompressible navier–stokes equations in time-varying complex geometries, Journal of Computational Physics 229 (12) (2010) 4476–4500.
  * [35] D.-V. Le, B. C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, Journal of Computational Physics 220 (1) (2006) 109–138.
  * [36] C. Brehm, C. Hader, H. F. Fasel, A locally stabilized immersed boundary method for the compressible navier–stokes equations, Journal of Computational Physics 295 (2015) 475–504.
  * [37] C. Liu, X. Zheng, C. Sung, Preconditioned multigrid methods for unsteady incompressible flows, Journal of Computational physics 139 (1) (1998) 35–57.
  * [38] C. H. Lee, A. J. Gil, G. Greto, S. Kulasegaram, J. Bonet, A new jameson–schmidt–turkel smooth particle hydrodynamics algorithm for large strain explicit fast dynamics, Computer Methods in Applied Mechanics and Engineering 311 (2016) 71–111.
  * [39] C. Zhang, Y. Zhu, Y. Yu, D. Wu, M. Rezavand, S. Shao, X. Hu, An artificial damping method for total lagrangian sph method with application in biomechanics, Engineering Analysis with Boundary Elements 143 (2022) 1–13.
  * [40] D. Wu, C. Zhang, X. Tang, X. Hu, An essentially non-hourglass formulation for total lagrangian smoothed particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 407 (2023) 115915.
  * [41] C. Zhang, J. Wang, M. Rezavand, D. Wu, X. Hu, An integrative smoothed particle hydrodynamics method for modeling cardiac function, Computer Methods in Applied Mechanics and Engineering 381 (2021) 113847.
  * [42] M. Aguirre, A. J. Gil, J. Bonet, A. A. Carreno, A vertex centred finite volume jameson–schmidt–turkel (jst) algorithm for a mixed conservation formulation in solid dynamics, Journal of Computational Physics 259 (2014) 672–699.