††thanks: Corresponding author<EMAIL_ADDRESS>
# Generalized and high-efficiency arbitrary-positioned buffer for smoothed
particle hydrodynamics

Shuoguo Zhang  Yu Fan  School of Engineering and Design, Technical University
of Munich, Garching, 85748, Germany  Yaru Ren  School of Engineering and
Design, Technical University of Munich, Garching, 85748, Germany  State Key
Laboratory of Hydraulics and Mountain River Engineering, Sichuan University,
Chengdu, Sichuan, China  Bin Qian  China Southwest Architecture Design and
Research Institute Corp. Ltd., Chengdu, Sichuan, China  Xiangyu Hu  School of
Engineering and Design, Technical University of Munich, Garching, 85748,
Germany

###### Abstract

This paper develops an arbitrary-positioned buffer for the smoothed particle
hydrodynamics (SPH) method, whose generality and high efficiency are achieved
through two techniques. First, with the local coordinate system established at
each arbitrary-positioned in-/outlet, particle positions in the global
coordinate system are transformed into those in it via coordinate
transformation. Since one local axis is located perpendicular to the
in-/outlet boundary, the position comparison between particles and the
threshold line or surface can be simplified to just this coordinate dimension.
Second, particle candidates subjected to position comparison at one specific
in-/outlet are restricted to those within the local cell-linked lists nearby
the defined buffer zone, which significantly enhances computational efficiency
due to a small portion of particles being checked. With this developed buffer,
particle generation and deletion at arbitrary-positioned in- and outlets of
complex flows can be handled efficiently and straightforwardly. We validate
the effectiveness and versatility of the developed buffer through 2-D and 3-D
non-/orthogonal uni-/bidirectional flows with arbitrary-positioned in- and
outlets, driven by either pressure or velocity boundary conditions.

††preprint: AIP/123-QED

##  I Introduction

When simulating engineering flows, flow paths typically do not repeat, and in-
and outlets may not align to facilitate the flow recirculation Holmes and
Pivonka (2021). For this kind of smoothed particle hydrodynamics (SPH) flow
simulation, the classic periodic boundary condition Braun  _et al._ (2015);
Negi and Ramachandran (2022a), where newly inflow particles are generated by
recycling outflow particles through translation, is incapable for particle
generation/deletion at in-/outlet. In comparison, the open boundary condition
Federico  _et al._ (2012); Vacondio  _et al._ (2012); Alvarado-Rodriguez  _et
al._ (2017); Ferrand  _et al._ (2017); Negi, Ramachandran, and Haftu (2020)
can theoretically handle this case easily due to its independent particle
handling (i.e., particle generation/deletion) at in-/outlet. Several methods
have been proposed in SPH literature regarding open boundary conditions. The
simplest method is to employ the buffer zone consisting of several layers of
buffer particles Braun  _et al._ (2015); Martin, Basa, and Quinlan (2009);
Tafuni  _et al._ (2018); Zhang  _et al._ (2023). At inlet, when a buffer
particle crosses the buffer-fluid interface, a duplicate particle with the
same physical state is generated, which then moves freely as an ordinary fluid
particle Wang  _et al._ (2019). Meanwhile, the original buffer particle is
recycled to a periodic position at the inlet boundary. At outlet, outflowing
particles are deleted once they cross the domain boundary Lyu  _et al._
(2022). A similar procedure, i.e., checking the particle’s relative position
to the prescribed threshold line (in two dimensions, 2-D) or surface (in three
dimensions, 3-D), is also employed in the mirror-axis technique Hirschler  _et
al._ (2016); Kunz  _et al._ (2016) for particle generation and deletion.
Monteleone et al. Monteleone, Monteforte, and Napoli (2017) utilized a conical
scan region with an angle $\beta$ for particle generation, comparing the angle
$\gamma$ between the cone axis and the line connecting the neighboring
particle with the cone vertex to $\beta$. If there are no effective particles
inside the conical scan region, i.e., $\gamma>\beta/2$ for all neighboring
particles, a new particle is created. However, when determining whether to
establish the conical scan region on a particle, the measured distance
$d_{ip}$ from this particle to the in-/outlet boundary needs to be compared
with both the particle spacing $dp$ and the smoothing length $h$, i.e.,
$dp<d_{ip}<h$. In contrast, the semi-analytical boundary treatment Ferrand
_et al._ (2017); Leroy  _et al._ (2016) does not rely on position comparison
for particle generation. Specifically, when the evolutive mass of fixed-
position inlet vertex particles reaches a certain threshold, a new fluid
particle is created at the same location. However, regarding particle
deletion, as same to the conical scan approach Monteleone, Monteforte, and
Napoli (2017), the semi-analytical boundary treatment still requires checking
whether outflowing particles are outside the domain through position
comparison.

Obviously, in these methods mentioned above, comparing the particle position
to a prescribed threshold line or surface is the basic and necessary operation
for particle generation and deletion. In SPH literature, the principles of
particle generation and deletion, along with most relevant test cases, are
typically illustrated based on simple orthogonal flows, where the in-/outflow
direction is parallel to one coordinate axis. As the in-/outlet boundary is
generally perpendicular to the in-/outflow direction, the prescribed threshold
line or surface for checking the particle’s relative position is usually
perpendicular to this coordinate axis as well. With this simplification, the
position comparison is performed in only one coordinate dimension. However, in
complex non-orthogonal flows, especially those with multiple arbitrary-
positioned in- and outlets, the above solution becomes non-trivial when
comparing the particle position to an arbitrary-positioned threshold line or
surface. Although directly determining the particle’s relative position
through the equation of the threshold line or surface can still work, the
complex procedure makes it inefficient and non-generalizable. In SPH
literature, some non-orthogonal flows with multiple arbitrary-positioned in-
and outlets have also been conducted Monteleone, Monteforte, and Napoli
(2017); Leroy  _et al._ (2016), but systematic elaborations on how to conduct
position comparison at arbitrary-positioned in-/outlet, especially with a
generalized and high-efficiency approach, have rarely been mentioned.
Furthermore, regarding particle detection, the straightforward way of
sequentially checking all fluid particle positions will significantly decrease
computational efficiency, especially in 3-D flows with multiple in- and
outlets.

In this paper, based on the buffer approach, we introduce the coordinate
transformation to develop a generalized and high-efficiency approach for
conducting position comparison in particle generation/deletion at arbitrary-
positioned in-/outlet. The particle positions in the global coordinate system
are first transformed into those in the local coordinate system established at
in-/outlet. Since one local axis is located perpendicular to the in-/outlet
boundary, their position comparison with the threshold line or surface is
performed only in this coordinate dimension. To further enhance computational
efficiency, at one specific in-/outlet, the particle candidates subjected to
position comparison are restricted to those within the local cell-linked lists
nearby the defined buffer zone. The remainder of this paper is organized as
follows. First, Section II provides a brief overview of the underlying weakly-
compressible SPH (WCSPH) method. In Section III, the rationale of arbitrary-
positioned uni-/bidirectional in-/outflow buffer is elaborated. Section IV
details pressure and velocity boundary conditions implemented at arbitrary-
positioned in-/outlet. The effectiveness and versatility of the developed
arbitrary-positioned buffer are demonstrated by several 2-D and 3-D flow
examples in Sections V and VI, including 2-D sloped-positioned VIPO (Velocity
Inlet, Pressure Outlet) and PIVO (Pressure Inlet, Velocity Outlet) channel
flows, 2-D T- and Y-shaped channel flows, the 2-D rotating sprinkler, the 3-D
sloped-positioned pulsatile pipe flow, and the 3-D aorta flow. Finally,
concluding remarks are given in Section VII. The code accompanying this work
is implemented in the open-source SPHinXsys library Zhang  _et al._ (2021a),
and is available at the project website https://www.sphinxsys.org and the
corresponding Github repository.

##  II Weakly-compressible SPH method

###  II.1 Governing equations

Within the Lagrangian framework, the mass- and momentum-conservation equations
are respectively written as

$\frac{d\rho}{dt}=-\rho\nabla\cdot\mathbf{v},$ (1)

and

$\frac{d\mathbf{v}}{dt}=-\frac{1}{\rho}\nabla p+\nu\nabla^{2}\mathbf{v},$ (2)

where $\rho$ is the density, $t$ the time, $\mathbf{v}$ the velocity, $p$ the
pressure, and $\nu$ the kinematic viscosity. Under the weakly-compressible
assumption, the system of Eqs. (1) and (2) are closed by the artificial
isothermal equation of state (EoS)

$p=c_{0}^{2}(\rho-\rho_{0}),$ (3)

where $c_{0}$ is the artificial sound speed, and $\rho_{0}$ the initial
reference density. In order to control the density variation around 1$\%$,
corresponding to the Mach number $M\approx 0.1$ suggested by Monaghan Monaghan
(1994), $c_{0}$ is set to be $10U_{max}$ with $U_{max}$ denoting the
anticipated maximum flow speed.

###  II.2 SPH discretization

Both the mass- and momentum-conservation equations Eqs. (1) and (2) are
discretized using the Riemann-based SPH scheme Zhang  _et al._ (2021a); Hu and
Adams (2006); Zhang, Hu, and Adams (2017a); Zhang, Rezavand, and Hu (2020);
Zhang  _et al._ (2021b), in respect to particle $i$

$\frac{d\rho_{i}}{dt}=2\rho_{i}\sum_{j}\frac{m_{j}}{\rho_{j}}(\mathbf{v}_{i}-\mathbf{v}^{*}_{ij})\cdot\nabla
W_{ij},$ (4)

and

$\frac{d\mathbf{v}_{i}}{dt}=-2\sum_{j}m_{j}(\dfrac{P^{*}_{ij}}{\rho_{i}\rho_{j}})\nabla
W_{ij}+2\sum_{j}m_{j}\frac{\eta\mathbf{v}_{ij}}{\rho_{i}\rho_{j}r_{ij}}\frac{\partial
W_{ij}}{\partial r_{ij}},$ (5)

where $m$ is the particle mass, $\eta$ the dynamic viscosity, subscript $j$
the neighbor particles, and $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$
the relative velocity between particles $i$ and $j$. Also, $\nabla W_{ij}$
denotes the gradient of the kernel function $W(|\mathbf{r}_{ij}|,h)$, where
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, and $h$ the smooth length.

Furthermore,
$\mathbf{v}^{*}=U^{*}\mathbf{e}_{ij}+(\mathbf{\overline{v}}_{ij}-\overline{U}\mathbf{e}_{ij})$
with $\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$ and
$\mathbf{\overline{v}}_{ij}=(\mathbf{v}_{i}+\mathbf{v}_{j})/2$ the average
velocity between particles $i$ and $j$. Here, constructed along the unit
vector $\mathbf{-e}_{ij}$ pointing from particle $i$ to $j$, the solutions
$U^{*}$ and $P^{*}$ of the inter-particle one-dimensional Riemann problem can
be given as Zhang, Hu, and Adams (2017b)

$\begin{cases}U^{*}=\overline{U}+\dfrac{P_{L}-P_{R}}{2\overline{\rho}c_{0}}\\\\[8.53581pt]
P^{*}=\overline{P}+\frac{1}{2}\overline{\rho}c_{0}(U_{L}-U_{R})\\\\[8.53581pt]
\end{cases},$ (6)

with the initial states of the Riemann problem

$\begin{cases}(\rho_{L},U_{L},P_{L})=(\rho_{i},\mathbf{-v}_{i}\cdot\mathbf{e}_{ij},p_{i})\\\\[8.53581pt]
(\rho_{R},U_{R},P_{R})=(\rho_{j},\mathbf{-v}_{j}\cdot\mathbf{e}_{ij},p_{j})\end{cases}.$
(7)

Here, $\overline{U}=(U_{L}+U_{R})/2$, $\overline{P}=(P_{L}+P_{R})/2$, and
$\overline{\rho}=(\rho_{L}+\rho_{R})/2$ are inter-particle averages. $L$ and
$R$ denotes the initial left and right states, respectively.

However, the utilization of the original intermediate pressure $P^{*}$ in Eq.
(6) during the SPH discretization may lead to an excessive dissipation. To
address this issue, a supplementary low dissipation Riemann solver Zhang, Hu,
and Adams (2017b) is employed, which incorporates a modification on $P^{*}$
while maintaining the intermediate velocity $U^{*}$ in Eq. (6) unconstrained

$P^{*}=\overline{P}+\frac{1}{2}\beta\overline{\rho}(U_{L}-U_{R}),$ (8)

where $\beta=min\>(3max(U_{L}-U_{R},0),c_{0})$ represents the low dissipation
limiter.

In order to increase the computational efficiency, the dual-criteria time
stepping method Zhang, Rezavand, and Hu (2020) is employed herein. In
addition, to decrease the accumulated density error during long-term
simulations, the particle density is re-initialized before being updated in
continuity equation Eq. (4) at each advection time step with Zhang  _et al._
(2021a)

$\rho_{i}=\rho_{0}\dfrac{\sum W_{ij}}{\sum W^{0}_{ij}},$ (9)

where the superscript $0$ represents the reference value in the initial
configuration.

###  II.3 Transport-velocity formulation

To avoid particle clumping and void regions in SPH simulations where the
tensile instability is present Lind  _et al._ (2012); Oger  _et al._ (2016),
regularization is implemented to maintain a uniform particle distribution
Litvinov, Hu, and Adams (2015); Negi and Ramachandran (2022b). Currently, the
particle shifting technique (PST) Lind  _et al._ (2012); Oger  _et al._
(2016); Negi and Ramachandran (2022b); Skillen  _et al._ (2013); Khayyer,
Gotoh, and Shimizu (2017); Sun  _et al._ (2017) and the transport-velocity
formulation (TVF) Zhang, Hu, and Adams (2017a); Adami, Hu, and Adams (2013)
are two typical schemes to address this issue. In the present work, the TVF
scheme is adopted, and the particle advection velocity
$\widetilde{\mathbf{v}}$ is rewritten as

$\widetilde{\mathbf{v}}_{i}(t+\delta t)=\mathbf{v}_{i}(t)+\delta
t\left(\frac{\widetilde{d}\mathbf{v}_{i}}{dt}-p_{max}\sum_{j}\frac{2m_{j}}{\rho_{i}\rho_{j}}\frac{\partial
W_{ij}}{\partial r_{ij}}\mathbf{e_{ij}}\right),$ (10)

where the global background pressure
$p_{max}=\lambda\rho_{0}\mathbf{v}_{max}^{2}$ with $\mathbf{v}_{max}$ denoting
the maximum particle velocity at each time step Zhang, Rezavand, and Hu
(2020). The empirical coefficient $\lambda$ is chosen within the range of 5.0
to 10.0 Zhang, Hu, and Adams (2017a).

Since the TVF scheme is applicable only to inner fluid particles far away from
in-/outlet boundaries Adami, Hu, and Adams (2013), accurately identifying
boundary particles is a necessary prior. To avoid misidentifying and missing
boundary particles Zhang  _et al._ (2023); Dilts (1999); Haque and Dilts
(2007); Lee  _et al._ (2008), the spatio-temporal identification approach
Zhang  _et al._ (2023) is utilized in present work, where approximately three
layers of fluid particles are identified as boundary particles.

##  III The rationale of arbitrary-positioned buffer

In this section, we introduce the rationale of arbitrary-positioned buffer
working for uni- and bidirectional flows, respectively. Although the schematic
provided here is based on channel flow, this arbitrary-positioned buffer can
also be implemented in free-surface open-channel flows with arbitrary-
positioned in- and outlets. In this work, the C2 Wendland kernel Wendland
(1995) with a compact support of $2h$ is employed. Since the present smoothing
length $h=1.3dp$ following the general practice, the buffer zone should
consist of at least three layers of particles to ensure full support for the
inner-domain particles next to the buffer.

###  III.1 Unidirectional in-/outflow buffer

Figure 1: The rationale of arbitrary-positioned unidirectional in-/outflow
buffer.

In Figure 1, a two-dimensional complex-geometric channel flow with multiple
arbitrary-positioned in- and outlets is sketched to demonstrate the rationale
of the arbitrary-positioned unidirectional in-/outflow buffer. Through the
translation
$\overrightarrow{OO^{{}^{\prime}}}=[X_{O^{{}^{\prime}}},Y_{O^{{}^{\prime}}}]^{T}$
and rotation $\theta$, the local coordinate system
$O^{{}^{\prime}}X^{{}^{\prime}}Y^{{}^{\prime}}$ is established at each
in-/outlet, with its axis direction
$\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}$ aligned to the in-/outflow
direction. Here, $[X_{O^{{}^{\prime}}},Y_{O^{{}^{\prime}}}]^{T}$ is the
position of origin $O^{{}^{\prime}}$ in the global coordinate system. The
buffer, with width $a$ and length $b$, is centered at the origin
$O^{{}^{\prime}}$ of the local coordinate system
$O^{{}^{\prime}}X^{{}^{\prime}}Y^{{}^{\prime}}$. As shown in Figure 1, one
side of the buffer is aligned with the in-/outlet boundary.

In the 2-D local coordinate system, the particle position
$\mathbf{r}^{{}^{\prime}}=(X^{{}^{\prime}},Y^{{}^{\prime}})^{T}$ could be
obtained through the following coordinate transformation

$\begin{bmatrix}X^{{}^{\prime}}\\\
Y^{{}^{\prime}}\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}\left(\begin{bmatrix}X\\\
Y\end{bmatrix}-\begin{bmatrix}X_{O^{{}^{\prime}}}\\\
Y_{O^{{}^{\prime}}}\end{bmatrix}\right),$ (11)

where $\mathbf{r}=[X,Y]^{T}$ is the corresponding particle position in the
global coordinate system, and the sign of $\theta$ depends on the
anticlockwise ($+$) or clockwise ($-$) rotation.

In comparison, for the coordinate transformation in the 3-D case, the unit
vector $\overrightarrow{\mathbf{n}}=\left[n_{x},n_{y},n_{z}\right]^{T}$ of the
rotation axis is first given by normalizing the cross product of axial vectors
$\overrightarrow{OX}$ and $\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}$ in
the global coordinate system, which is a vector perpendicular to the plane
spanned by $\overrightarrow{OX}$ and
$\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}$ with the direction given by
the right-hand rule

$\overrightarrow{\mathbf{n}}=\dfrac{\overrightarrow{OX}\times\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}}{\left|\overrightarrow{OX}\times\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}\right|},$
(12)

and then its antisymmetric matrix $\mathbf{S}$ can be written as

$\mathbf{S}=\begin{bmatrix}0&-n_{z}&n_{y}\\\ n_{z}&0&-n_{x}\\\
-n_{y}&n_{x}&0\end{bmatrix}.$ (13)

According to the Rodrigues’ rotation formula Cheng and Gupta (1989); Fraiture
(2009); Dai (2015), the rotation matrix $\mathbf{R}$ is defined as

$\mathbf{R}=\mathbf{I}+\mathbf{S}\sin\theta+\mathbf{S}^{2}(1-\cos\theta),$
(14)

where $\mathbf{I}$ is 3$\times$3 identity matrix, and the sign of $\theta$ is
positive.

Consequently, the particle position
$\mathbf{r}^{{}^{\prime}}=(X^{{}^{\prime}},Y^{{}^{\prime}},Z^{{}^{\prime}})^{T}$
in the 3-D local coordinate system could be given as

$\begin{bmatrix}X^{{}^{\prime}}\\\ Y^{{}^{\prime}}\\\
Z^{{}^{\prime}}\end{bmatrix}=\mathbf{R}\left(\begin{bmatrix}X\\\ Y\\\
Z\end{bmatrix}-\begin{bmatrix}X_{O^{{}^{\prime}}}\\\ Y_{O^{{}^{\prime}}}\\\
Z_{O^{{}^{\prime}}}\end{bmatrix}\right).$ (15)

In the 2-D inflow buffer, buffer particles will be periodically moved back to
the inlet bound with
$\mathbf{r}^{{}^{\prime}}=\mathbf{r}^{{}^{\prime}}-\left[a,0\right]^{T}$
Tafuni  _et al._ (2018); Zhang  _et al._ (2023), once they cross the buffer-
fluid interface. Here,
$\mathbf{r}^{{}^{\prime}}=\mathbf{r}^{{}^{\prime}}-\left[a,0,0\right]^{T}$ in
the 3-D case. Consequently, due to its fixed buffer particle configuration,
buffer particles should be identified from all fluid particles in the
computational domain but only for once. In detail, by checking
$\mathbf{r}^{{}^{\prime}}$ of all fluid particles before the flow simulation,
they will be identified as buffer particles if both
$\left|X^{{}^{\prime}}\right|<a/2$ and $\left|Y^{{}^{\prime}}\right|<b/2$ (an
additional criterion of $\left|Z^{{}^{\prime}}\right|<c/2$ for 3-D case, with
$c$ denoting the buffer height), and then added into the set of buffer
particles. Note that, each inflow buffer has its own set of buffer particles.
During the subsequent flow evolution, when these identified buffer particles
cross the buffer-fluid interface with $X^{{}^{\prime}}>a/2$ (same for 3-D
case), duplicate particles are generated with the same states as the crossing
ones, and then treated as ordinary fluid particles in the inner flow domain.
Meanwhile, the original buffer particles will be recycled to periodic
positions at inlet bound, and their new positions in the global coordinate
system are obtained through the inverse transformation of coordinates

$\begin{bmatrix}X\\\ Y\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}^{T}\left(\begin{bmatrix}X^{{}^{\prime}}\\\
Y^{{}^{\prime}}\end{bmatrix}-\begin{bmatrix}a\\\
0\end{bmatrix}\right)+\begin{bmatrix}X_{O^{{}^{\prime}}}\\\
Y_{O^{{}^{\prime}}}\end{bmatrix}.$ (16)

Correspondingly, for the 3-D case, the above inverse transformation of
coordinates is modified as

$\begin{bmatrix}X\\\ Y\\\
Z\end{bmatrix}=\mathbf{R}^{T}\left(\begin{bmatrix}X^{{}^{\prime}}\\\
Y^{{}^{\prime}}\\\ Z^{{}^{\prime}}\end{bmatrix}-\begin{bmatrix}a\\\ 0\\\
0\end{bmatrix}\right)+\begin{bmatrix}X_{O^{{}^{\prime}}}\\\
Y_{O^{{}^{\prime}}}\\\ Z_{O^{{}^{\prime}}}\end{bmatrix}.$ (17)

Algorithm 1 summarizes the working procedure of the unidirectional inflow
buffer described above.

Data: $i$: all fluid particles in the computational domain, $j$: buffer
particles, $n$: the number of particles being checked.

1

2Procedure Buffer particle identification $\hfill\triangleright$Execute only
once

3 for  _$i=1$ to $n$_ do

4 if  _$\left|X_{i}^{{}^{\prime}}\right| <\frac{a}{2}$ &
$\left|Y_{i}^{{}^{\prime}}\right|<\frac{b}{2}$ for 2-D case &
$\left|Z_{i}^{{}^{\prime}}\right|<\frac{c}{2}$ for 3-D case_ then

5 add it into the set of buffer particles

6 end if

7

8 end for

9Procedure Particle generation

10 for  _$j=1$ to $n$_ do

11 if  _$X_{j}^{{}^{\prime}} >\frac{a}{2}$ for both 2-D and 3-d cases_ then

12 duplicate it, and then recycle it through the corresponding invese
transformation of coordinates Eq. 16 or Eq. 17

13 end if

14

15 end for

Algorithm 1 Working procedure of the arbitrary-positioned unidirectional
inflow buffer.

For the 2-D outflow buffer, as shown in Figure 1, only the fluid particles
within the nearby local cell-linked lists at each time step are checked. Note
that, the coverage area of the local cell-linked lists is a little larger than
the size of outflow buffer, approximately one layer of cells. If
$X^{{}^{\prime}}>a/2$ of these checked particles (same for 3-D case), they
will be directly deleted as outflow ones. Furthermore, due to the newly inflow
fluid particles and outflow buffer particles, the buffer particle
configuration of the outflow buffer necessitates an update at each time step.
In detail, through the relabeling approach Zhang  _et al._ (2024), at the end
of each time step (after the cell-linked lists have been rebuilt), if both
$\left|X^{{}^{\prime}}\right|<a/2$ and $\left|Y^{{}^{\prime}}\right|<b/2$ (an
additional criterion of $\left|Z^{{}^{\prime}}\right|<c/2$ for 3-D case) of
these checked particles, they will be relabeled as buffer particles by
assigning the same integer, otherwise, they will not. Note that, each outflow
buffer can have its exclusive integer for relabeling. Algorithm 2 summarizes
the working procedure of the unidirectional outflow buffer described above.

Data: $i$: fluid particles within the local cell-linked lists, $n$: the number
of particles being checked.

1

2Procedure Particle deletion

3 for  _$i=1$ to $n$_ do

4 if  _$X_{i}^{{}^{\prime}} >\frac{a}{2}$ for both 2-D and 3-d cases_ then

5 delete it

6 end if

7

8 end for

9Procedure Buffer particle identification (after the cell-linked lists have
been rebuilt at the end of present time step)

10 for  _$i=1$ to $n$_ do

11 if  _$\left|X_{i}^{{}^{\prime}}\right| <\frac{a}{2}$ &
$\left|Y_{i}^{{}^{\prime}}\right|<\frac{b}{2}$ for 2-D case &
$\left|Z_{i}^{{}^{\prime}}\right|<\frac{c}{2}$ for 3-D case_ then

12 relabel it as buffer particle

13 end if

14

15 end for

Algorithm 2 Working procedure of the arbitrary-positioned unidirectional
outflow buffer.

###  III.2 Bidirectional in-/outflow buffer

Figure 2: The rationale of arbitrary-positioned bidirectional buffer.

Based on Figure 1, Figure 2 comparatively illustrates the different setting
and rationale of the arbitrary-positioned bidirectional buffer in the
bidirectional flow. Here, the axis direction
$\overrightarrow{O^{{}^{\prime}}X^{{}^{\prime}}}$ of the established local
coordinate system at each in-/outlet is always aligned to the inflow
direction, as shown in Figure 2.

Similarly, due to the newly inflow fluid particles and outflow buffer
particles, the buffer particle configuration of the bidirectional buffer
should also be updated at the end of each time step in the same way of the
arbitrary-positioned unidirectional outflow buffer, where the particles being
checked are also restricted to those within the nearby local cell-linked
lists. Differently, particle generation and deletion during each time step
only focus on identified buffer particles. If their $X^{{}^{\prime}}>a/2$
(same for 3-D case), the same procedure for particle generation as in the
arbitrary-positioned unidirectional inflow buffer will be executed.
Conversely, if their $X^{{}^{\prime}}<-a/2$ (same for 3-D case), they will be
directly deleted. Algorithm 3 summarizes the working procedure of the
arbitrary-positioned bidirectional buffer described above.

Data: $i$: fluid particles within the local cell-linked lists, $j$: buffer
particles, $n$: the number of particles being checked.

1 Procedure Particle generation and deletion

2 for  _$j=1$ to $n$_ do

3 if  _$X_{j}^{{}^{\prime}} >\frac{a}{2}$ for both 2-D and 3-d cases_ then

4 duplicate it, and then recycle it through the corresponding invese
transformation of coordinates Eq. 16 or Eq. 17

5 end if

6

7 if  _$X_{j}^{{}^{\prime}} <-\frac{a}{2}$ for both 2-D and 3-d cases_ then

8 delete it

9 end if

10

11 end for

12Procedure Buffer particle identification (after the cell-linked lists have
been rebuilt at the end of present time step)

13 for  _$i=1$ to $n$_ do

14 if  _$\left|X_{i}^{{}^{\prime}}\right| <\frac{a}{2}$ &
$\left|Y_{i}^{{}^{\prime}}\right|<\frac{b}{2}$ for 2-D case &
$\left|Z_{i}^{{}^{\prime}}\right|<\frac{c}{2}$ for 3-D case_ then

15 relabel it as buffer particle

16 end if

17

18 end for

Algorithm 3 Working procedure of the arbitrary-positioned bidirectional
buffer.

##  IV In-/outlet boundary conditions

In this section, we introduce how to implement pressure and velocity boundary
conditions at arbitrary-positioned in-/outlet through buffer particles. To
cope with the truncation for density summation, the density reinitialization
of Eq. (9) is not carried out for near-boundary particles (i.e., buffer
particles herein).

###  IV.1 Pressure boundary condition

In the present work, the pressure boundary condition follows the work of Zhang
et al. Zhang  _et al._ (2024). In the global coordinate system, the
discretized momentum equation Eq. (5), when applied to buffer particles, is
modified with the prescribed boundary pressure $p_{b}$ incorporated

$\begin{split}\frac{d\mathbf{v}_{i}}{dt}=&-2\sum_{j}m_{j}(\dfrac{P^{*}_{ij}}{\rho_{i}\rho_{j}})\nabla
W_{ij}+2p_{b}\sum_{j}(\dfrac{m_{j}}{\rho_{i}\rho_{j}})\nabla W_{ij}\\\
&+2\sum_{j}m_{j}\frac{\eta\mathbf{v}_{ij}}{\rho_{i}\rho_{j}r_{ij}}\frac{\partial
W_{ij}}{\partial r_{ij}}\end{split}.$ (18)

As mentioned by Holmes et al. Holmes and Pivonka (2021), the velocity updated
by Eq. (18) should be further projected to the normal direction $\hat{u}$ of
in-/outlet boundary in the global coordinate system

$\mathbf{v}_{i}=(\mathbf{v}_{i}\cdot\hat{u})\hat{u}.$ (19)

Here, $\hat{u}$ could be expressed in the 2-D case as

$\hat{u}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}^{T}\begin{bmatrix}1\\\
0\end{bmatrix}=\begin{bmatrix}\cos\theta\\\ \sin\theta\end{bmatrix},$ (20)

and in the 3-D case as

$\hat{u}=\mathbf{R}^{T}\begin{bmatrix}1\\\ 0\\\ 0\end{bmatrix}.$ (21)

Furthermore, to satisfy the consistency condition at the boundary, the density
of the newly populated (actually recycled) particles in the buffer is obtained
following the boundary pressure and EoS as

$\rho_{i}=\rho_{0}+p_{b}/c_{0}^{2}.$ (22)

###  IV.2 Velocity boundary condition

Figure 3: The coordinate transformation on pre-defined velocity profile.

For the sake of simplicity, the pre-defined velocity profile is first
implemented on buffer particles within the local coordinate system, for
example in Figure 3,

$\begin{bmatrix}\mathbf{v}_{X^{{}^{\prime}}}\\\
\mathbf{v}_{Y^{{}^{\prime}}}\end{bmatrix}=\begin{bmatrix}1-(\frac{Y^{{}^{\prime}}}{R})^{2}\\\
0\end{bmatrix}\qquad\text{in 2-D},$ (23)

and

$\begin{bmatrix}\mathbf{v}_{X^{{}^{\prime}}}\\\
\mathbf{v}_{Y^{{}^{\prime}}}\\\
\mathbf{v}_{Z^{{}^{\prime}}}\end{bmatrix}=\begin{bmatrix}1-(\frac{Y^{{}^{\prime}}}{R})^{2}\\\
0\\\ 0\end{bmatrix}\qquad\text{in 3-D},$ (24)

where $R$ is the channel radius.

Then the corresponding velocity distribution in the global coordinate system
can be derived through coordinate transformation on buffer particle velocity

$\begin{bmatrix}\mathbf{v}_{X}\\\
\mathbf{v}_{Y}\end{bmatrix}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}^{T}\begin{bmatrix}\mathbf{v}_{X^{{}^{\prime}}}\\\
\mathbf{v}_{Y^{{}^{\prime}}}\end{bmatrix}\qquad\text{in 2-D},$ (25)

and

$\begin{bmatrix}\mathbf{v}_{X}\\\ \mathbf{v}_{Y}\\\
\mathbf{v}_{Z}\end{bmatrix}=\mathbf{R}^{T}\begin{bmatrix}\mathbf{v}_{X^{{}^{\prime}}}\\\
\mathbf{v}_{Y^{{}^{\prime}}}\\\
\mathbf{v}_{Z^{{}^{\prime}}}\end{bmatrix}\qquad\text{in 3-D}.$ (26)

The pressure boundary condition by Eq.(18) should also be imposed at the
boundary to eliminate the truncated error in approximating pressure gradient,
but the prescribed boundary pressure $p_{b}$ in Eq.(18) is given as $p_{i}$,
i.e., $p_{b}=p_{i}$. Furthermore, both density and pressure of newly populated
(actually recycled) particles in the buffer remain unchanged.

##  V 2-D Numerical tests

In this section, to validate the effectiveness of the developed arbitrary-
positioned buffer in the 2-D case, several 2-D flows are simulated, including
sloped-positioned VIPO and PIVO channel flows, T- and Y-shaped channel flows,
and rotating sprinkler. For the treatment of wall boundary, a one-sided
Riemann solver is employed herein Zhang, Hu, and Adams (2017b).

###  V.1 Sloped-positioned VIPO and PIVO channel flows

Figure 4: Schematics of sloped-positioned VIPO (top panel) and PIVO (bottom
panel) channel flows. Velocity profile is extracted at the cross-section
$\mathbf{A-A}$ in the simulation.

As shown in Figure 4, 2-D sloped-positioned VIPO (Velocity Inlet, Pressure
Outlet) and PIVO (Pressure Inlet, Velocity Outlet) channel flows with the
available analytical solution are first simulated herein. At the velocity
in-/outlet, the prescribed velocity profile
$\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}})$ in the local coordinate system
is analytically determined by Holmes and Pivonka (2021); Morris, Fox, and Zhu
(1997); Takeda, Miyama, and Sekiya (1994); Sigalotti  _et al._ (2003); Holmes,
Williams, and Tilke (2011)

$\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}})=\frac{\mathtt{\Delta}P}{2\eta
L}(\frac{d}{2}-y^{{}^{\prime}})(\frac{d}{2}+y^{{}^{\prime}}),$ (27)

where $y^{{}^{\prime}}\in[-\frac{d}{2},\frac{d}{2}]$ with $d=0.001m$ denoting
the distance between plates, and $\mathtt{\Delta}P=0.1Pa$ the pressure drop
across a flow length of $L=0.004m$. The dynamic viscosity is set as
$\eta=\sqrt{\rho d^{3}\mathtt{\Delta}P/8LRe}$ with the fluid density
$\rho=1000kg/m^{3}$ and the Reynolds number $Re=50$. The fluid and solid-wall
particles are initialized on the Cartesian Lattice with a uniform particle
spacing of $dp=d/50$.

Figure 5: Pressure contours of VIPO (top panel) and PIVO (bottom panel)
channel flows. Time instant $t=\infty$.

At the pressure outlet of VIPO flow and the pressure inlet of PIVO flow, the
prescribed boundary pressure $p_{b}$ is set to $0Pa$ and $0.1Pa$,
respectively. Consequently, in the VIPO flow, the predicted boundary pressure
at the velocity inlet should be $0.1\text{ Pa}$, while at the velocity outlet
of PIVO flow, it should be kept at $0\text{ Pa}$. Figure 5 gives the pressure
contours of VIPO and VIPO flows at time $t=\infty$. It can be obviously
observed that the pressure distributions align well with the above-mentioned
theoretical expectations.

Figure 6: Comparison of numerical and analytical velocity profiles in the
local coordinate system for VIPO (left panel) and PIVO (right panel) channel
flows. Time instant $t=\infty$. The RMSEP errors for VIPO and PIVO channel
flows are $1.91\%$ and $0.96\%$, respectively.

Figure 6 illustrates the comparison between numerical and analytical velocity
profiles for VIPO and PIVO channel flows. The distribution of velocity data
points consistently exhibits a parabolic shape, closely aligning with the
analytical solution. Specifically, the maximum velocity occurs at the
channel’s centerline, while the velocity diminishes to zero near the walls due
to the no-slip boundary condition enforced on the solid walls Adami, Hu, and
Adams (2012). The numerical accuracy is further measured with the Root Mean
Square Error Percentage (RMSEP) over the whole velocity profile Holmes and
Pivonka (2021)

$RMSEP=\sqrt{\dfrac{1}{N}\sum_{n=1}^{N}\left(\dfrac{\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})-\mathbf{\widetilde{v}}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})}{\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})}\right)^{2}},$
(28)

where $N$ is the number of measuring points,
$\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})$ and
$\mathbf{\widetilde{v}}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})$ the analytical
and numerical velocities, respectively. The RMSEP errors of VIPO and PIVO
channel flows are $1.91\%$ and $0.96\%$, respectively, demonstrating a high
level of concordance with the analytical solution.

###  V.2 T- and Y-shaped channel flows

Here we conduct two preliminary 2-D application cases where daughter pipes are
connected to the main pipe at different angles, i.e., T- and Y-shaped channel
flows in Figure 7, and compare them to the results obtained from Fluent
simulations. In the local coordinate system, the inflow velocity profile is
given as $\mathbf{v}_{x^{{}^{\prime}}}=400(0.05^{2}-(y^{{}^{\prime}})^{2})$
with the characteristic velocity $U_{max}=1m/s$. Hence, the dynamic viscosity
$\eta=0.1\rho U_{max}/Re$ with the Reynolds number $Re=100$ and fluid density
$\rho=1000kg/m^{3}$. The boundary pressure $p_{b}$ at the pressure outlet is
set to $0pa$. The uniform particle/grid resolution $dp$ is $10^{-3}m$.

Figure 7: T- and Y-shaped channel flow schematics with the highlighted
unidirectional inflow buffer. Velocity profile is extracted at the cross-
section $\mathbf{A-A}$ in the simulation. Figure 8: Qualitative comparison of
velocity contours from SPH and Fluent simulations for T- (top panel) and Y-
(bottom panel) shaped channel flows. Time instant $t=\infty$.

The left panels of Figure 8 illustrate the SPH velocity contours of T- and
Y-shaped channel flows at time $t=\infty$, respectively. At the fork junction,
i.e., the stagnation point, the fluid velocity decreases to almost zero,
especially in the T-shaped channel flow. The right panels of Figure 8 show the
corresponding results of Fluent simulation with the same model settings.
Although good agreement is already demonstrated in the qualitative comparison
of Figure 8, we further conduct a quantitative comparison with velocity
profiles extracted at the outlet, as shown in Figure 9. Both the shapes of
velocity profiles from SPH and Fluent simulations present consistent parabolic
forms, with nearly identical velocity magnitudes. By replacing the analytical
velocity $\mathbf{v}_{x^{{}^{\prime}}}(y^{{}^{\prime}}_{n})$ in Eq. (28) as
the numerical velocity from the Fluent simulation, the RMSEP errors of T- and
Y-shaped channel flows are measured as $2.29\%$ and $3.01\%$, respectively,
which are considered acceptable.

Figure 9: Comparison of SPH and Fluent velocity profiles in the local
coordinate system for T- (left panel) and Y- (right panel) shaped channel
flows. Time instant $t=\infty$. $R$ denotes the radius of daughter channel.
The RMSEP errors for T- and Y-shaped channel flows are $2.29\%$ and $3\%$,
respectively.

###  V.3 Rotating sprinkler

To better demonstrate the versatility of the developed arbitrary-positioned
buffer in engineering applications, the 2-D rotating sprinkler, referring to a
lawn sprinkler, is simulated here, which is realized by rotating buffer
particles with a constant angular velocity $\omega$, as depicted in Figure 10.
In this top view simulation, gravity $\mathbf{g}$ is ignored, and the fluid is
considered inviscid with the density of $\rho=1000kg/m^{3}$. The inflow
velocity profile in the sprinkler is defined within the local coordinate
system as $\mathbf{v}_{x^{{}^{\prime}}}=2500(0.02^{2}-(y^{{}^{\prime}})^{2})$.
The uniform particle spacing is $dp=0.005m$.

Figure 10: Top view of rotating sprinkler. The sprinkler consisting of buffer
particles is rotating with the constant angular velocity $\omega$. Figure 11:
The flow chart on rotating the sprinkler through coordinate transformation.

Figure 11 illustrates the flowchart for rotating the sprinkler through
coordinate transformation. The rotation angle $\theta$ is given as
$\theta=\omega t+\theta_{0}$, where $\theta_{0}$ denotes the initial rotation
angle at time $t=0$. In the present case, the angular velocity $\omega=2rad/s$
and $\theta_{0}=0$. Before rotating the sprinkler at $t_{1}$, the particle
generation in this unidirectional inflow buffer has already been completed in
the established local coordinate system. Then, based on this local coordinate
system, the positions and velocities of buffer particles are transformed to
new ones through coordinate rotation at $t_{1}$. After this rotation, the
local coordinate system is reestablished according to the new position of
sprinkler, preparing it for particle generation in next time step. Similarly,
based on this newly established local coordinate system, the sprinkler
undergoes the next rotation at $t_{2}$, following the same procedure as at
$t_{1}$. Thus, the iterative cycle continues.

Figure 12: Velocity contours of rotating sprinkler at different time instants.

Figure 12 shows the velocity contours at different time instants. During the
continuous rotation, the velocity profile within the sprinkler consistently
maintains the prescribed parabolic shape, and new inflow particles at the
sprinkler’s outlet are effectively generated. Due to inertia, the fluid bulk
exhibits a curved shape, consistent with physical principles. At these four
instants, all rotated angles of the sprinkler match the predicted angles based
on the angular velocity. All these flow behaviors demonstrate the good
reliability and performance of the arbitrary-positioned buffer developed by
introducing the coordinate transformation.

##  VI 3-D Numerical tests

In this section, to validate the effectiveness of the developed arbitrary-
positioned buffer in the 3-D case, two 3-D flows with the same wall boundary
condition as the 2-D flow are simulated, including the sloped-positioned
pulsatile pipe flow and aorta flow.

###  VI.1 Sloped-positioned pulsatile pipe flow

Building upon the 2-D sloped-positioned VIPO channel flow in section V.1, an
additional dimension, typically representing the radial direction in the 3-D
case, has been included to model the 3-D sloped-positioned pulsatile pipe
flow, while all other physical and geometrical parameters remain unchanged.
Differently, the pressure boundary condition is imposed at both in- and
outlets. When the constant pressure gradient is modified to follow a cosine
function over time, i.e.,
$P^{{}^{\prime}}=\mathtt{\Delta}P/L=0.1\cos(t)/0.004=25\cos(t)$, as
illustrated in the upper panel of Figure 13, a bidirectional pulsatile flow is
facilitated. The analytical velocity profile in the local coordinate system,
i.e., $\mathbf{v}_{x^{{}^{\prime}}}(r,t)$ at time $t$ and radius $r$ from the
pipe axis, could be given as Womersley (1955)

$\mathbf{v}_{x^{{}^{\prime}}}(r,t)=\mathbf{Re}\bigg{\\{}\sum_{n=0}^{N}\dfrac{iP_{n}^{{}^{\prime}}}{\rho
n\omega}\bigg{[}1-\dfrac{J_{0}(\alpha
n^{\frac{1}{2}}i^{\frac{3}{2}}\frac{r}{R})}{J_{0}(\alpha
n^{\frac{1}{2}}i^{\frac{3}{2}})}\bigg{]}e^{in\omega t}\bigg{\\}},$ (29)

where $r\in(0,R)$, $\alpha=R\sqrt{\omega\rho/\eta}$ the dimensionless
Womersley number, $\omega$ the angular frequency of the first harmonic of a
Fourier series of an oscillatory pressure gradient, $P_{n}^{{}^{\prime}}$ the
pressure gradient magnitude for the frequency $n\omega$, $J_{0}(\cdot)$ the
Bessel function of first kind and order zero, $i$ the imaginary number, and
$\mathbf{Re}\left\\{\cdot\right\\}$ the real part of a complex number.

In the present pulsatile pipe flow, the Womersley number $\alpha=0.8409$ is
less than 2, which determines a parabolic shape of the pulsatile flow profile.
Moreover, during the time evolution, consistent with the cosine-varying
pressure gradient, both velocity magnitude and direction exhibit periodicity.
The lower panel of Figure 13 correspondingly shows this periodic evolution of
the parabolic shape by extracting velocity profiles at 13 instants between
$t=2\pi$ and $4\pi$. During the initial quarter, the parabolic velocity
profile maintains a forward flow direction, but with a diminishing magnitude
over time. The flow undergoes a reversal at $t=2.5\pi$, followed by gradual
acceleration until reaching the maximum velocity at $t=3\pi$. In the latter
half of the period, this flow behavior repeats, but with the opposite
evolution of flow direction. Moreover, employing RMSEP for accuracy assessment
reveals errors of approximately $3\%$ across all 13 time instants,
demonstrating a robust agreement with their corresponding analytical
solutions.

Furthermore, aligning with the analytical prediction of mixed in- and outflows
at time $t=11.116s$ shown in the bottom panel of Figure 14, the top panel of
Figure 14 depicts the mixed forward and reverse flows occurring at the
in-/outlet boundary of the practical pulsatile pipe flow. This demonstrates
the good performance of the arbitrary-positioned bidirectional buffer in
managing particle generation and deletion.

Figure 13: The pressure gradient as a cosine function (upper panel) in the
$-45^{\circ}$ sloped-positioned pulsatile pipe flow, and the comparison (lower
panel) of numerical and analytical velocity profiles in the local coordinate
system at 13 indicative instants. The horizontal grid dimension in the lower
panel spans a velocity range from $0$ to 5.0E-3 $m/s$. The errors at 13
instants, from $t=2\pi$ to $4\pi$, are $1.43\%$, $2.27\%$, $2.8\%$, $3.45\%$,
$3.29\%$, $2.38\%$, $2.91\%$, $2.65\%$, $2.89\%$, $3.55\%$, $2.43\%$, $1.85\%$
and $1.35\%$, respectively.

Figure 14: Velocity profile with mixed forward and reverse flows at the
in-/outlet boundary (top panel) of $-45^{\circ}$ sloped-positioned pulsatile
pipe flow, corresponding to the analytical velocity profile (bottom panel) in
the local coordinate system at time instant $t=11.116s$. The velocity vector
in top panel is scaled to a uniform length, with the flow direction indicated
by the arrow.

###  VI.2 Aorta flow

To investigate the feasibility of the developed arbitrary-positioned buffer
for handling practical 3-D complex flows with multiple arbitrary-positioned
in- and outlets, a preliminary flow simulation is conducted on an aorta model,
as shown in Figure 15. The vascular wall is treated as a rigid body, while the
blood is characterized as a Newtonian fluid with a density of 1060 $kg/m^{3}$
Laganà  _et al._ (2005) and a viscosity of 0.00355 $Pa\cdot s$ Gallo  _et al._
(2012).

Figure 15: Schematic of the aorta flow.

At the inlet, a velocity boundary condition featuring a plug velocity profile
is enforced to facilitate unidirectional inflow. However, the velocity
magnitude exhibits pulsatility, as depicted in the top panel of Figure 17,
which can be approximated as Fourier series in term of time Taib  _et al._
(2019)

$\mathbf{v}_{inlet}(t)=a_{0}+\sum_{n=1}^{8}a_{n}cos(nt\omega)+\sum_{n=1}^{8}b_{n}sin(nt\omega),$
(30)

where the empirical variables for the Fourier series are listed in Table 1.

Table 1: The empirical variables for the Fourier series Taib  _et al._ (2019) Parameter | Value | Parameter | Value
---|---|---|---
$a_{0}$ | 0.3782 | $\omega$ | 8.302
$a_{1}$ | -0.1812 | $b_{1}$ | -0.07725
$a_{2}$ | 0.1276 | $b_{2}$ | 0.01466
$a_{3}$ | -0.08981 | $b_{3}$ | 0.04295
$a_{4}$ | 0.04347 | $b_{4}$ | -0.06679
$a_{5}$ | -0.05412 | $b_{5}$ | 0.05679
$a_{6}$ | 0.02642 | $b_{6}$ | -0.01878
$a_{7}$ | 0.008946 | $b_{7}$ | 0.01869
$a_{8}$ | -0.009005 | $b_{7}$ | -0.01888

At all outlets, the pressure boundary condition is applied by employing the
time-dependent boundary pressure $P(t)$, thus the bidirectional flow happens
at outlets during the cyclical blood flow, as shown in bottom panel of Figure
16. Here, the 3-element (RCR) Windkessel model de Almeida Jansen Catanho,
Sinha, and Vijayan (2012) is employed for the real-time computing of $P(t)$

$(1+\dfrac{R_{p}}{R_{d}})Q(t)+CR_{p}\frac{dQ(t)}{dt}=\dfrac{P(t)}{R_{d}}+C\dfrac{dP(t)}{dt},$
(31)

where $Q(t)$ is the blood flow volume, $C$ the arterial compliance, $R_{p}$
the characteristic resistance, and $R_{d}$ the peripheral resistance. Table 2
lists the values of Windkessel parameters Kim (2007); Sudharsan and Cherry
Kemmerling (2018) for allocation across diverse daughter vessels.

Table 2: Parameters for the Windkessel outlet boundary conditions Kim (2007); Sudharsan and Cherry Kemmerling (2018) Vessel | $R_{p}(dynes\>s/cm^{5})$ | $C(cm^{5}/dynes)$ | $R_{d}(dynes\>s/cm^{5})$
---|---|---|---
Right common carotid | 1180 | 7.70E-5 | 18400
Right subclavian | 1040 | 8.74E-5 | 16300
Left common carotid | 1180 | 7.70E-5 | 18400
Left subclavian | 970 | 9.34E-5 | 15200
Descending aorta | 188 | 4.82E-4 | 2950

The top and middle panels of Figure 16 shows the velocity contour and slices
during peak systole. In the present plug-inlet boundary condition, the
velocity peak isn’t always at the center of the cross-section. Especially in
the descending aorta, due to inertia, the flow inside the curved vessel is
pushed towards the outer side of the arch, which is in good agreement with the
observations in literature Sudharsan and Cherry Kemmerling (2018).
Furthermore, as depicted in the bottom panel of Figure 16, which indicates the
corresponding velocity vectors of boundary particles, the present arbitrary-
positioned bidirectional buffer can effectively handle particle generation and
deletion at these typical mixed in- and outflow boundaries.

Figure 16: Velocity contour and slices of aorta flow during the peak systole,
along with indicating corresponding velocity vectors of boundary particles.
Time instant $t=4.9s$. The velocity vector in bottom panel is scaled to a
uniform length, with the flow direction indicated by the arrow.

Figure 17: The time history of the inflow velocity and all outlet relative
pressure.

Figure 17 illustrates the temporal evolution of relative pressure at all
outlets across 7 cardiac cycles. Notably, all observed outlet relative
pressure waveforms exhibit a characteristic folding line shape, consistent
with established analytical and numerical findings in literature de Almeida
Jansen Catanho, Sinha, and Vijayan (2012). It is pertinent to acknowledge
that, relative to real physiological conditions, the parameters governing
inflow velocity, Windkessel parameters of all daughter vessels, and the
geometric specifications of the aorta stem from disparate data sources,
thereby lacking strict alignment. Moreover, the physiological pressure data
corresponding to the present aorta model are also unavailable. Consequently,
our primary emphasis lies in assessing the fidelity of the blood pressure
waveform rather than pressure magnitudes.

##  VII Conclusion

In this paper, we develop a generalized and high-efficiency arbitrary-
positioned uni-/bidirectional buffer for SPH method. By employing coordinate
transformation, the position comparison for particle generation and deletion
is conducted in the local coordinate system established at each in-/outlet.
Since one local axis is located perpendicular to the in-/outlet boundary, the
position comparison can be simplified and conducted only in this coordinate
dimension. To further enhance computational efficiency, the particle
candidates subjected to position comparison are restricted to those within the
local cell-linked lists nearby the defined buffer zone. The numerical results
of 2-D and 3-D flows demonstrate the excellent effectiveness and versatility
of the developed buffer for particle generation and deletion at arbitrary-
positioned in- and outlets. As the future work, we will implement the
developed buffer in various free-surface open-channel flows with multiple
arbitrary-positioned in- and outlets.

## References

  * Holmes and Pivonka (2021) D. W. Holmes and P. Pivonka, “Novel pressure inlet and outlet boundary conditions for smoothed particle hydrodynamics, applied to real problems in porous media flow,” Journal of Computational Physics 429, 110029 (2021).
  * Braun  _et al._ (2015) S. Braun, L. Wieth, R. Koch, and H. Bauer, “A framework for permeable boundary conditions in SPH: Inlet, Outlet, Periodicity,” (2015).
  * Negi and Ramachandran (2022a) P. Negi and P. Ramachandran, “How to train your solver: Verification of boundary conditions for smoothed particle hydrodynamics,” Physics of Fluids 34 (2022a).
  * Federico  _et al._ (2012) I. Federico, S. Marrone, A. Colagrossi, F. Aristodemo, and M. Antuono, “Simulating 2d open-channel flows through an sph model,” European Journal of Mechanics - B/Fluids 34, 35–46 (2012).
  * Vacondio  _et al._ (2012) R. Vacondio, B. Rogers, P. Stansby, and P. Mignosa, “SPH modeling of shallow flow with open boundaries for practical flood simulation,” Journal of Hydraulic Engineering 138, 530–541 (2012).
  * Alvarado-Rodriguez  _et al._ (2017) C. Alvarado-Rodriguez, J. Klapp, L. Sigalotti, J. Dominguez, and E. Sanchez, “Nonreflecting outlet boundary conditions for incompressible flows using SPH,” Computers and Fluids 159, 177–188 (2017).
  * Ferrand  _et al._ (2017) M. Ferrand, A. Joly, C. Kassiotis, D. Violeau, A. Leroy, F. Morel, and B. Rogers, “Unsteady open boundaries for SPH using semi-analytical conditions and riemann solver in 2D,” Computer Physics Communications 210, 29–44 (2017).
  * Negi, Ramachandran, and Haftu (2020) P. Negi, P. Ramachandran, and A. Haftu, “An improved non-reflecting outlet boundary condition for weakly-compressible SPH,” Computer Methods in Applied Mechanics and Engineering 367, 113119 (2020).
  * Martin, Basa, and Quinlan (2009) L. Martin, M. Basa, and N. Quinlan, “Permeable and non-reflecting boundary conditions in SPH,” International Journal for Numerical Methods in Fluids 61, 709–724 (2009).
  * Tafuni  _et al._ (2018) A. Tafuni, J. Domínguez, R. Vacondio, and A. Crespo, “A versatile algorithm for the treatment of open boundary conditions in smoothed particle hydrodynamics GPU models,” Computer Methods in Applied Mechanics and Engineering 342, 604–624 (2018).
  * Zhang  _et al._ (2023) S. Zhang, W. Zhang, C. Zhang, and X. Hu, “A lagrangian free-stream boundary condition for weakly compressible smoothed particle hydrodynamics,” Journal of Computational Physics 490, 112303 (2023).
  * Wang  _et al._ (2019) P. Wang, A. Zhang, F. Ming, P. Sun, and H. Cheng, “A novel nonreflecting boundary condition for fluid dynamics solved by smoothed particle hydrodynamics,” Journal of Fluid Mechanics 860, 81–114 (2019).
  * Lyu  _et al._ (2022) H.-G. Lyu, P. Sun, A. Colagrossi, and A.-M. Zhang, “Towards sph simulations of cavitating flows with an eosb cavitation model,” Acta Mechanica Sinica 39 (2022).
  * Hirschler  _et al._ (2016) M. Hirschler, P. Kunz, M. Huber, F. Hahn, and U. Nieken, “Open boundary conditions for isph and their application to micro-flow,” Journal of Computational Physics 307, 614–633 (2016).
  * Kunz  _et al._ (2016) P. Kunz, M. Hirschler, M. Huber, and U. Nieken, “Inflow/outflow with dirichlet boundary conditions for pressure in isph,” Journal of Computational Physics 326, 171–187 (2016).
  * Monteleone, Monteforte, and Napoli (2017) A. Monteleone, M. Monteforte, and E. Napoli, “Inflow/outflow pressure boundary conditions for smoothed particle hydrodynamics simulations of incompressible flows,” Computers and Fluids 159, 9–22 (2017).
  * Leroy  _et al._ (2016) A. Leroy, D. Violeau, M. Ferrand, L. Fratter, and A. Joly, “A new open boundary formulation for incompressible sph,” Computers & Mathematics with Applications 72, 2417–2432 (2016).
  * Zhang  _et al._ (2021a) C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, J. Wang, and X. Hu, “SPHinXsys: An open source multi-physics and multi-resolution library based on smoothed particle hydrodynamics,” Computer Physics Communications 267, 108066 (2021a).
  * Monaghan (1994) J. Monaghan, “Simulating free surface flows with SPH,” Journal of Computational Physics 110, 399–406 (1994).
  * Hu and Adams (2006) X. Hu and N. Adams, “A multi phase SPH method for macroscopic and mesoscopic flows,” Journal of Computational Physics 213, 844–861 (2006).
  * Zhang, Hu, and Adams (2017a) C. Zhang, X. Hu, and N. Adams, “A generalized transport-velocity formulation for smoothed particle hydrodynamics,” Journal of Computational Physics 337, 216–232 (2017a).
  * Zhang, Rezavand, and Hu (2020) C. Zhang, M. Rezavand, and X. Hu, “Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics,” Journal of Computational Physics 404, 109135 (2020).
  * Zhang  _et al._ (2021b) C. Zhang, Y. Wei, F. Dias, and X. Hu, “An efficient fully lagrangian solver for modeling wave interaction with oscillating wave surge converter,” Ocean Engineering 236, 109540 (2021b).
  * Zhang, Hu, and Adams (2017b) C. Zhang, X. Hu, and N. Adams, “A weakly compressible SPH method based on a low-dissipation Riemann solver,” Journal of Computational Physics 335 (2017b).
  * Lind  _et al._ (2012) S. Lind, R. Xu, P. Stansby, and B. Rogers, “Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves.” Journal of Computational Physics 231, 1499–1523 (2012).
  * Oger  _et al._ (2016) G. Oger, S. Marrone, D. Le Touzé, and M. de Leffe, “Sph accuracy improvement through the combination of a quasi-lagrangian shifting transport velocity and consistent ale formalisms,” Journal of Computational Physics 313, 76–98 (2016).
  * Litvinov, Hu, and Adams (2015) S. Litvinov, X. Hu, and N. A. Adams, “Towards consistence and convergence of conservative sph approximations,” Journal of Computational Physics 301, 394–401 (2015).
  * Negi and Ramachandran (2022b) P. Negi and P. Ramachandran, “Techniques for second-order convergent weakly compressible smoothed particle hydrodynamics schemes without boundaries,” Physics of Fluids 34, 087125 (2022b).
  * Skillen  _et al._ (2013) A. Skillen, S. Lind, P. Stansby, and B. Rogers, “Incompressible smoothed particle hydrodynamics (sph) with reduced temporal noise and generalised fickian smoothing applied to body-water slam and efficient wave-body interaction,” Computer Methods in Applied Mechanics and Engineering 265, 163–173 (2013).
  * Khayyer, Gotoh, and Shimizu (2017) A. Khayyer, H. Gotoh, and Y. Shimizu, “Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in isph context,” Journal of Computational Physics 332, 236–256 (2017).
  * Sun  _et al._ (2017) P. Sun, A. Colagrossi, S. Marrone, and A. Zhang, “The delta plus-sph model: Simple procedures for a further improvement of the sph scheme,” Computer Methods in Applied Mechanics and Engineering 315, 25–49 (2017).
  * Adami, Hu, and Adams (2013) S. Adami, X. Hu, and N. Adams, “A transport-velocity formulation for smoothed particle hydrodynamics,” Journal of Computational Physics 241, 292–307 (2013).
  * Dilts (1999) G. Dilts, “Moving-least-squares-particle hydrodynamics—I. Consistency and stability,” International Journal for Numerical Methods in Engineering 44, 1115–1155 (1999).
  * Haque and Dilts (2007) A. Haque and G. Dilts, “Three-dimensional boundary detection for particle methods,” Journal of Computational Physics 226, 1710–1730 (2007).
  * Lee  _et al._ (2008) E. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, and P. Stansby, “Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method,” Journal of Computational Physics 227, 8417–8436 (2008).
  * Wendland (1995) H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree,” Advances in Computational Mathematics 4, 389–396 (1995).
  * Cheng and Gupta (1989) H. Cheng and K. C. Gupta, “An Historical Note on Finite Rotations,” Journal of Applied Mechanics 56, 139–145 (1989).
  * Fraiture (2009) L. Fraiture, “A history of the description of the three-dimensional finite rotation,” The Journal of the Astronautical Sciences 57, 207–232 (2009).
  * Dai (2015) J. S. Dai, “Euler–rodrigues formula variations, quaternion conjugation and intrinsic connections,” Mechanism and Machine Theory 92, 144–152 (2015).
  * Zhang  _et al._ (2024) S. Zhang, Y. Fan, D. Wu, C. Zhang, and X. Hu, “Dynamical pressure boundary condition for weakly-compressible smoothed particle hydrodynamics,” arXiv preprint arXiv:2403.09485 (2024).
  * Morris, Fox, and Zhu (1997) J. P. Morris, P. J. Fox, and Y. Zhu, “Modeling Low Reynolds Number Incompressible Flows Using SPH,” Journal of Computational Physics 136, 214–226 (1997).
  * Takeda, Miyama, and Sekiya (1994) H. Takeda, S. M. Miyama, and M. Sekiya, “Numerical Simulation of Viscous Flow by Smoothed Particle Hydrodynamics,” Progress of Theoretical Physics 92, 939–960 (1994).
  * Sigalotti  _et al._ (2003) L. D. G. Sigalotti, J. Klapp, E. Sira, Y. Meleán, and A. Hasmy, “SPH simulations of time-dependent Poiseuille flow at low Reynolds numbers,” Journal of Computational Physics 191, 622–638 (2003).
  * Holmes, Williams, and Tilke (2011) D. Holmes, J. Williams, and P. Tilke, “Smooth particle hydrodynamics simulations of low reynolds number flows through porous media,” International Journal for Numerical and Analytical Methods in Geomechanics 35, 419 – 437 (2011).
  * Adami, Hu, and Adams (2012) S. Adami, X. Hu, and N. Adams, “A generalized wall boundary condition for smoothed particle hydrodynamics,” Journal of Computational Physics 231, 7057–7075 (2012).
  * Womersley (1955) J. R. Womersley, “Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known,” The Journal of Physiology 127, 553–563 (1955).
  * Laganà  _et al._ (2005) K. Laganà, R. Balossino, F. Migliavacca, G. Pennati, E. L. Bove, M. R. de Leval, and G. Dubini, “Multiscale modeling of the cardiovascular system: application to the study of pulmonary and coronary perfusions in the univentricular circulation,” Journal of Biomechanics 38, 1129–1141 (2005).
  * Gallo  _et al._ (2012) D. Gallo, G. De Santis, F. Negri, D. Tresoldi, R. Ponzini, D. Massai, M. Deriu, P. Segers, B. Verhegghe, G. Rizzo, and U. Morbiducci, “On the use of in vivo measured flow rates as boundary conditions for image-based hemodynamic models of the human aorta: Implications for indicators of abnormal flow,” Annals of biomedical engineering 40, 729–41 (2012).
  * Taib  _et al._ (2019) I. Taib, M. S. A. Paisal, A. E. Ismail, A. M. Tajul Arifin, and N. Darlis, “An analysis of blood pressure waveform using windkessel model for normotensive and hypertensive conditions in carotid artery,” Journal of Advanced Research in Fluid Mechanics and Thermal Sciences 57, 69–85 (2019).
  * de Almeida Jansen Catanho, Sinha, and Vijayan (2012) M. T. de Almeida Jansen Catanho, M. Sinha, and V. Vijayan, “Model of aortic blood flow using the windkessel effect,” (2012).
  * Kim (2007) _On coupling a lumped-parameter heart model with a three-dimensional finite element model of the aorta_ , ASME Summer Bioengineering Conference (2007).
  * Sudharsan and Cherry Kemmerling (2018) M. Sudharsan and E. Cherry Kemmerling, “The effect of inlet and outlet boundary conditions in image-based cfd modeling of aortic flow,” BioMedical Engineering OnLine 17 (2018).