Dual-criteria time stepping for weakly compressible smoothed particle hydrodynamics

TL;DR

This paper introduces a dual-criteria time-stepping method for weakly compressible smoothed particle hydrodynamics that enhances computational efficiency by decoupling particle interaction updates from pressure relaxation steps.

Contribution

It proposes a novel dual-criteria time-stepping approach that improves efficiency in WCSPH simulations by using separate advection and acoustic time criteria.

Findings

Reduced CPU cost in simulations

Maintained accuracy in benchmark tests

Potential for fluid-structure interaction modeling

Abstract

Implementing particle-interaction configuration and time integration are performance intensive essentials of particle-based methods. In this paper, a dual-criteria time-stepping method is proposed to improve the computational efficiency of the weakly-compressible smoothed particle hydrodynamic (WCSPH) method for modeling incompressible flows. The key idea is to introduce an advection time criterion, which is based on fluid velocity field, for recreating the particle-interaction configuration. Within this time criterion, several steps of pressure relaxation determined by the acoustic time criterion, based on the artificial speed of sound, can be carried out without updating the particle interaction configuration and much larger time-step sizes compared with the conventional counterpart. The method has shown optimized computational performance through CPU cost analysis. Good accuracy and performance is obtained for the presented benchmarks implying promising potential of the proposed method for incompressible flow and fluid-structure interaction simulations.

Tables

Table 1. Table 1: CPU time for simulating the dam-break flow with three spatial resolutions.

	Total CPU time (s)
Number of Particles	64560 (H/dp = 20)	269472 (H/dp = 40)	1290364 (H/dp = 80)
Traditonal	1683.89	29217.4	525621
Present	583.45	10408.2	187314
Acceleration ratio	2.88	2.81	2.80

Table 2. Table 2: Flow-induced vibration of a beam attached to a cylinder: comparison of the present computational results against previous studies.

Refs.	Amplitude in $y$ direction	Frequency
Turek and Hron [49]	$0.83$	$0.19$
Bhardwaj and Mittal [50]	$0.92$	$0.19$
Tian et al. [51]	$0.784$	$0.19$
Present method	$0.93$	$0.178$

Equations

$$\{\begin{array}{l}\frac{\text{d}\rho }{\text{d}t}=-\rho \nabla \cdot \mathbf{v}\\ \rho \frac{\text{d}\mathbf{v}}{\text{d}t}=-\nabla p+\eta {\nabla }^2\mathbf{v}+\rho \mathbf{g}\end{array},$$

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

$$\frac{\text{d}{\rho }_i}{\text{d}t}={\rho }_i\sum _j\frac{m_j}{{\rho }_j}{v}_{ij}\cdot {\nabla }_i{W}_{ij}=2{\rho }_i\sum _j\frac{m_j}{{\rho }_j}(v_i-{\overline{v}}_{ij})\cdot {\nabla }_i{W}_{ij}.$$

$$\frac{\text{d}\mathbb{v}_{i}}{\text{d}t}=-\sum_{j}m_{j}\bigg{(}\frac{p_{i}+p_{j}}{\rho_{i}\rho_{j}}\bigg{)}\nabla_{i}W_{ij}=-2\sum_{j}m_{j}\frac{\overline{p}_{ij}}{\rho_{i}\rho_{j}}\nabla_{i}W_{ij},$$

$$\begin{split}\bigg{(}\frac{\text{d}\mathbb{v}_{i}}{\text{d}t}\bigg{)}^{(\nu)}&=-\sum_{j}m_{j}\alpha\frac{hc}{\overline{\rho}}\frac{\mathbb{v}_{ij}\cdot\mathbb{r}_{ij}}{|\mathbb{r}_{ij}|^{2}}\nabla_{i}W_{ij}\end{split}$$

$${\overline{v}}_{ij}=U^{\ast }{e}_{ij}+({\overline{v}}_{ij}-\overline{U}{e}_{ij}),{\overline{p}}_{ij}=P^{\ast },$$

$$\bigg{(}\frac{\text{d}\mathbb{v}_{i}}{\text{d}t}\bigg{)}^{(\nu)}=2\sum_{j}m_{j}\frac{\eta}{\rho_{i}\rho_{j}}\frac{{{\mathbf{v}}_{ij}}}{{{r}_{ij}}}\frac{\partial{{W}_{ij}}}{\partial{{r}_{ij}}},$$

$$\frac{\text{d}{\tilde{v}}_i}{\text{d}t}=\frac{\text{d}{v}_i}{\text{d}t}-2\sum _j{m}_j\frac{p^0}{{\rho }_i{\rho }_j}{\nabla }_i{W}_{ij},$$

$$\Delta t_{ad}={CFL}_{ad}\min \left(\frac{h}{|\mathbf{v}{|}_{max}},\frac{h^2}{\nu }\right),$$

$$\Delta t_{ac}={CFL}_{ac}\frac{h}{c+|\mathbf{v}{|}_{max}},$$

$${\rho }_i=\max ({\rho }^{\ast },{\rho }^0\frac{\sum W_{ij}}{\sum W_{ij}^0}),$$

$${\rho }_i={\rho }^0\frac{\sum W_{ij}}{\sum W_{ij}^0}.$$

$$\mathbf{v}_{i}^{n+\frac{1}{2}}=\mathbf{v}_{i}^{n}+\frac{1}{2}\Delta t_{ac}\big{(}\frac{d\mathbf{v}_{i}}{dt}\big{)}^{n}.$$

$$\begin{cases}\mathbf{r}_{i}^{n+1}=\mathbf{r}_{i}^{n+\frac{1}{2}}+\Delta t_{ac}\mathbf{v}_{i}^{n+\frac{1}{2}}\\ \rho_{i}^{n+1}=\rho_{i}^{n}+\frac{1}{2}\Delta t_{ac}\big{(}\frac{d\rho_{i}}{dt}\big{)}^{n+\frac{1}{2}}.\end{cases}$$

$$\mathbf{v}_{i}^{n+1}=\mathbf{v}_{i}^{n}+\frac{1}{2}\Delta t_{ac}\big{(}\frac{d\mathbf{v}_{i}}{dt}\big{)}^{n+1}.$$

$$P_s=\frac{{\sum }_s{p}_f{W}_{sf}}{{\sum }_s{W}_{sf}+ϵ},$$

$${\mathbf{v}}_i\leftarrow {\mathbf{v}}_i\alpha +\mathbf{v}({\mathbf{r}}_i)(1-\alpha ),$$

Full text

1 Introduction

As a fully Lagrangain particle-based method, smoothed particle hydrodynamics (SPH) was originally proposed by Lucy [1] and Gingold and Monaghan [2] for modeling astrophysics problems. Beside its original application [3, 4], SPH has been successfully exploited in a broad variety of problems ranging from solid mechanics [5, 6, 7, 8] to fluid mechanics [9, 10, 11] and fluid-structure interactions (FSI) [12, 13, 14]. Due to its Lagrangian feature, SPH is particularly well suited for modeling flows involving significantly varying topology and free surface [15, 16, 17, 18, 19].

In the context of modeling incompressible flows, the weakly-compressible SPH (WCSPH) method is widely applied and many studies have been devoted to portray the applicability of WCSPH to simulate scientific and engineering problems. Typical examples include wave breaking or overtopping [20, 21], wave impact [22, 13], violent sloshing [23, 24], as well as FSI problems [12]. A well-known limitation of WCSPH is the computational efficiency, in particular arising when one aims at long time simulations of large scale systems within reasonable computing cost. One type of techniques to alleviate this limitation is applying High Performance Computing (HPC) techniques by using Central Processing Units (CPUs), Graphics Processing Units (GPUs) and Many Integrated Cores (MICs) systems to accelerate the simulations [25, 26, 27]. Another type of the acceleration techniques is to modify the numerical algorithms to circumvent or optimize the critical computationally expensive operations [28, 29, 30].

SPH algorithms are based on pairwise interactions between neighboring particles through a Gaussian-like kernel function, which is radial-symmetric and has compact support. Therefore, implementing the particle-interaction configuration, i.e. determining the neighboring particles (by neighbor search) and computing corresponding kernel function values, is a critical aspect of the high-performance SPH solvers. Two different approaches, namely, cell-linked list (CLL) [28] and Verlet list (VL) [31], are widely used in the SPH community for implementing the particle-interaction configuration.

The CLL approach first partitions the computational domain into equisized cells, typically with the size of the support or cut-off radius of the kernel function, such that a CLL is created for each cell with all particles within the cell. Having the CLLs created, neighbor search is restricted to the nearest neighboring cells, 9 and 27 cells in two and three dimensions, respectively. Then, together with neighbor-searching operations, the CLL approach computes the kernel function values and interaction forces with the particles belong to these neighboring cells only if they are found within the cut-off radius. Note that since the CLL approach does not store the neighboring-particle identities and the corresponding kernel function values, usually the neighbor search is carried out multiple times, together with the kernel function values recomputed accordingly, during a single time step. For example, when the kick-drift-kick time integration algorithm is applied in a typical WCSPH method [32], neighbor search and computing the kernel values are carried out twice for one single time-step.

The VL approach also first divides the computational domain into equisized cells but typically with the size larger than the cut-off radius. For each particle, all particles within the adjacent cells are then checked to create a VL containing the references to all potential neighboring particles. In the simplest VL approach, since the cell size is the same as the cut-off radius, the VL contains only the real neighbors with non-vanishing kernel function values [29]. Note that a VL may be used for multiple times without executing neighbor search if the particle-interaction configuration can be obtained without updating the VL. Again, for example, in the above mentioned WCSPH method [32], when the simplest VL approach is applied, a VL can be used twice for one single time step.

However, as pointed by Dominguez et al. [29], the simplest VL approach does not show notable improvement on computational efficiency as the VL is only reused once. In order to increase the reuse of VL, the VL approach has to use cells with considerably larger size than the cut-off radius. Dominguez et al. [29] showed that with a 50% increase of the cell size, a performance gain of $8\%$ is achieved when the VL is reused for $13$ times in $7$ time steps. The reason for only incremental performance gain is that all the potential neighboring particles need to be checked for computing interaction force even only a small portion of them are real neighbors. This deficiency is due to the amount of potential particles scales to the cell size with the power of dimensions.

In this paper, we present a dual-criteria time stepping for the WCSPH method to improve computational efficiency when the simplest VL approach is employed. The basic idea is, apart from the acoustic time criterion as in previous WCSPH methods, to introduce an extra advection time criterion which is similar to that used in incompressible SPH (ISPH) method [33]. While the simplest VL is employed, such choice of time criteria leads to a dual-loop time-stepping algorithm in which the recreation of VLs and computation of kernel function values are carried out in the outer advection time-step loops and the pressure relaxation is carried out in the inner acoustic time-step loops. Note that, in comparison to ISPH, the present method circumvents the solution of pressure Poisson equation (PPE), which requires global quantities for iterations and converged PPE solution to ensure numerical stability. The reason for improved computational efficiency by the present method is twofold. First, larger time-steps compared with previous counterparts are permissible for the acoustic criterion. Second, even the VL contains only the real neighbors, it can be reused for many time during the acoustic time steps. The remainder of this paper is arranged as follows: Section 2 details the proposed time-steeping method. A number of validation test cases are presented and discussed together with computational efficiency analysis in Section 3 and concluding remarks of the present study are summarized in Section 4.

2 Method

In the Lagrangian frame, the conservation of mass and momentum can be written as

[TABLE]

where $\mathbf{v}$ is the velocity, $\rho$ the density, $p$ the pressure, $\eta$ the dynamic viscosity, $\mathbf{g}$ the gravity and $\frac{\text{d}}{\text{d}t}=\frac{\partial}{\partial t}+\mathbf{v}\cdot\nabla$ stands for material derivative. Following the weakly-compressible assumption [9, 34] for modeling incompressible flow, Eq. (1) is closed by an artificial isothermal equation of state (EoS)

[TABLE]

Eq. 2 assumes that the density varies around $1\%$ [34] if an artificial sound speed of $c=10U_{max}$ is employed, with $U_{max}$ being the maximum anticipated velocity within the flow.

2.1 WCSPH method

Following Refs. [11, 35, 36, 37], the WCSPH discretization of the continuity equation reads

[TABLE]

Here, $m_{j}$ is the particle mass, and $\mathbb{v}_{ij}=\mathbb{v}_{i}-\mathbb{v}_{j}$ and $\overline{\mathbb{v}}_{ij}=(\mathbb{v}_{i}+\mathbb{v}_{j})/2$ are the relative and average velocities between particle $i$ and $j$, respectively. $\nabla_{i}W_{ij}$ represents the gradient of the kernel function $W(|\mathbb{r}_{ij}|,h)$, where $\mathbb{r}_{ij}=\mathbb{r}_{i}-\mathbb{r}_{j}$ and $h$ is the smoothing length, with respect to particle $i$. Similarly, the discretization of the inviscid momentum equation gives

[TABLE]

where $\overline{p}_{ij}=(p_{i}+p_{j})/2$ denotes the average pressure between particle $i$ and $j$. To stabilize the inviscid simulation, an artificial viscosity term [36, 38] can be included as

[TABLE]

where $\overline{\rho}=(\rho_{i}+\rho_{j})/2$, and $0\leq\alpha\leq 1.0$ is a tunable parameter. An alternative approach to stabilize the simulation is introducing a Riemann solver to realize the pairwise particle interactions [37, 32]. Specifically, the inter-particle averages $\overline{\mathbb{v}}_{ij}$ and $\overline{p}_{ij}$ in Eqs. (3) and (4) are replaced by

[TABLE]

where $\overline{U}=\mathbb{v}_{ij}\cdot\mathbb{e}_{ij}$, $U^{\ast}$ and $P^{\ast}$ are the solution of inter-particle Riemann problem along the unit vector $\mathbb{e}_{ij}=-\mathbb{r}_{ij}/r_{ij}$ pointing form particle $i$ to particle $j$. For viscous flows, Eq. (5) is replaced by the physical shear force term obtained through [11]

[TABLE]

where $\eta$ is the dynamic viscosity of the fluid. The WCSPH method may suffer from particle clumping and void regions in high Reynolds number flows. As a remedy, the transport-velocity formulation [39, 17] with the following form

[TABLE]

is employable for modeling internal flows to improve the accuracy and stability. Here, $p^{0}$ is the background pressure and $\mathbb{\widetilde{v}}$ represents the particle transport velocity,

2.2 Dual-criteria time stepping

The dual-criteria time-stepping employs two time-step size criteria characterized by the particle advection and the acoustic velocities, respectively. The time-step size determined by the advection criterion, termed $\Delta t_{ad}$, has the following form

[TABLE]

where $CFL_{ad}=0.25$, $|\mathbf{v}|_{max}$ is the maximum particle advection velocity in the flow and $\nu$ the kinematic viscosity. Note that this criterion is the same as that used in the ISPH method [40]. The time-step size according to the acoustic criterion, termed $\Delta t_{ac}$, has the form

[TABLE]

where ${CFL}_{ac}=0.6$. Note that this criterion gives much larger time-step size than that employed in conventional time integration for WCSPH simulations [9, 35, 36].

In this paper, while the advection criterion controls the updating frequency of the simplest VL and the corresponding kernel function values, the acoustic criterion determines the frequency of the pressure relaxation process, i.e. the time integration of the particle density, pressure and velocity. Accordingly, during one advection step, the pressure relaxation process is carried out approximately $k\simeq\frac{\Delta t_{ad}}{\Delta t_{ac}}$ times. During the pressure relaxation process, the particle-interaction configuration is considered to be fixed in space, same as the ISPH method and temporally similar to the case for a Eulerian method. As a consequence, a large $CFL_{ac}=0.6$ number typically for a Eulerian method, for example the WENO scheme [41], is allowable without introducing numerical instability.

The details of the time stepping procedure are given in the following. At the beginning of the advection step, the fluid density field of free-surface flows is reinitialized by

[TABLE]

where $\rho^{*}$ denotes the density before re-initialization and superscript [math] represents the initial reference value. For flows without free surface, Eq. 11 is merely modified as

[TABLE]

Similar to the previous re-initialization approaches in Refs. [10, 17], Eq. 11 or 12 stabilizes the density which is updated by Eq. 3 in the pressure relaxation process without updating the particle-interaction configuration. Also, the viscous force is computed here and the transport-velocity formulation is applied if necessary. Having the time-step size $\Delta t_{ad}$ calculated, the pressure relaxation process is repeated using a velocity-Verlet scheme [31, 42] with the time-step size of $\Delta t_{ac}$, until the accumulated time interval is larger than $\Delta t_{ad}$. Here, we denote the beginning of an acoustic time step by superscript $n$, at the mid-point by $n+\frac{1}{2}$ and eventually at the new time-step by $n+1$. In the Verlet scheme, the velocity field is first updated to the mid-point value by

[TABLE]

Then particle position and density are updated to the next time step in the following form

[TABLE]

Finally, the velocity field is updated at the end of the time step as

[TABLE]

An overview of the proposed time-stepping method is shown in Fig. 1. Note that the present time stepping recovers the traditional scheme by removing the outer loop and density reinitialization, applying ${CFL}_{ac}=0.25$ and updating the VL and the kernel function values at every time step.

3 Validation tests

In this section, to validate the present method, a set of benchmark cases are simulated and the results are compared with those of previous experiments and simulations. These cases include several three-dimensional free-surface flows exhibiting violent phenomena such as wave impacts and breaking, and a two-dimensional FSI problem involving an elastic structure oscillating in flow. The $5th$-order Wendland kernel [43] with a smoothing length of $h=1.3dp$, where $dp$ is the initial particle spacing, and a cut-off radius of $2.6dp$ is employed in all the following simulations. As for the treatment of wall boundary, a one-sided Riemann solver is employed and we refer to Ref. [32] for more details.

3.1 Dam-break flow

In the first case, we consider a dam-break flow which is widely used as a benchmark and has experimental data [44] available for quantitative comparison. Following Lobovsky et al. [44], the schematic and setup parameters are given in Fig. 2. We consider an inviscid flow with a density of $\rho^{0}=1000~{}kg/m^{3}$ and a gravitational constant of $g=9.8~{}m/s^{2}$. According to the shallow water theory, the maximum velocity is estimated as $v_{max}=2\sqrt{gH}$ to set up the speed of sound, where $H$ is the initial water depth. Note that, for simplicity, the water column starts with zero initial pressure for the computation instead of being released by lifting up a holding gate quickly as in the experimental setup.

Fig. 3 illustrates several snapshots of the time evolution of free surface. Similar to previous numerical results [45, 32], the main features, viz. high roll-up along the downstream wall, a large reflected jet and a free surface breaking due to the re-entry of the backward wave, are well captured by the present method.

The predicted surge-wave front with increasing spatial resolutions is given in Fig. 4. It can be observed that the surge-wave front converges asymptotically to the shallow-water solution at later time. Note that, the analytical prediction is not applicable to the initial time, as the assumption of shallow-water equation does not hold. Similar to several previous simulations [42, 45, 13, 32], it is noted that numerically predicted surge-wave propagation is slightly faster than the experimental data, which is attributed to several uncertainties, e.g. turbulence, effect of air, wall roughness and repeatability of the experiment.

In the simulation, the measured pressure is obtained by kernel averaging from the fluid particles within the support radius, e.g, $2h$, of a specific sensor with

[TABLE]

where $\epsilon=1.0\times 10^{-6}$ is a small number to avoid zero denominator. Comparison of the time history of the predicted pressure signals against experimental data is plotted in Fig. 5. Considering that latter is post-processed by averaging the values of pressure on a sensor with a diameter of $d=4.2~{}mm$, a good agreement is noted except for the fluctuations, which are decreased with increasing the spatial resolution. Note that, the predicted pressure signal at $P3$ has lower magnitude than the experiment, but matches the occurrence time $t(\sqrt{g/H})\sim 2.8$ quite well. Similar discrepancies on pressure magnitude have been observed in previous numerical results of Cercos [47] and Zhang et al. [32], simulated by different WCSPH methods.

Here, we compare the computational efficiency of the present method with the traditional counterpart as interpreted in Section 2.2. The computations are carried out on an Intel(R) Xeon(R) L5520 2.27GHz Desktop computer with 24GiB RAM and Scientific Linux system (6.9). We choose the dam-break flow here because its has been used for performance evaluation in several previous studies [30, 47, 29]. To analyze the performance in strong scaling behavior, we evaluate the CPU wall-clock time for sequential computations of three different spatial resolutions until the dimensionless time of $10$. As shown in Tab. 1, an approximated speedup of $2.80$ is achieved by the present time-stepping method. This result indicates that, although the VL and corresponding kernel function values need more memory than the conventional CLL, it is a good choice for large scale simulations when such memory requirement is not a burden.

3.2 Liquid sloshing

Following Rafiee et al. [24], we consider an energetic violent free-surface flow, which involves liquid sloshing inside a rectangular tank with a low filling level. The geometric and setup parameters of the problem are given in Fig. 6. The tank motion is defined by a sinusoidal excitation in $x$-direction of $x=A\sin(2.0f\pi t)$, where $A=0.1~{}m$ and $f=0.496~{}s^{-1}$ are the amplitude and frequency, respectively. As $f$ is close to the natural frequency of the liquid, thereby resulting in a highly nonlinear violent free-surface motion. Initially, the particles are placed on a regular lattice with a particle spacing of $dp=L/260$.

Fig. 7 shows several snapshots of particle and pressure distributions. The main features, e.g., a traveling wave with a crest, the resulting bore breaks and impacts at the wall and then forming a high run-up, are well captured by the present method. Fig. 8 shows the zoom-in view of the free surface breaks and impacting during which a void region is generated and collapsed. This phenomenon is also observed in the experiments of Rafiee et al. [24], where the void region is filled with air.

The time history of the predicted impact pressure in comparison with experimental data is given in Fig. 9. A good agreement is noted, however, some oscillations are observed during the decay time and the impact pressure is overestimated. The discrepancies may be attributed to the neglected air phase in the present simulation. In the experimental investigations [24], there exists entrapped air pockets when the wave impacting occurs as shown in Fig. 8. The air pockets might cushion the impact and influence the peak pressure. Note that considerably more oscillatory impact pressure is obtained in the results of Rafiee et al. [24] than that predicted by the present simulation.

3.3 Dam-break flow with an obstacle

We consider a dam-break flow impacting a cuboid obstacle located on the downstream horizontal bed. This test case was first experimentally studied by Kleefsman et al. [48] with comparison against numerical results obtained by a Eulerian volume-of-fluid (VOF) method. This benchmark case has been frequently used to validate the applicability of particle-based methods in modeling FSI problems where the structure elasticity is neglected [32, 13]. The geometric and setup parameters are briefly described in Fig. 10. We consider an inviscid flow with a density of $\rho^{0}=1000~{}kg/m^{3}$ and the gravitational constant of $g=9.8~{}m/s^{2}$. The particles are initially placed on a regular lattice with a particle spacing of $dp=H/55$.

Fig. 11 shows several snapshots of the particle and pressure distributions at different time instants. It is observed that the free-surface profiles herein are in good agreement with the experimental and numerical results presented in [48], as well as those from previous simulations [32, 13]. Note that, similar to the experimental observations [48], a splash-up is produced after the surge wave impacts at the obstacle. Fig. 12 plots the time history of the pressure measured at $P1$ and $P2$, which indicate a good agreement with the experimental data except for the overestimated peak magnitude and the oscillations during pressure decay. These discrepancies are deemed to be reasonable due to the weakly-compressible formulation and the inviscid model employed in the present simulation.

3.4 Flow-induced vibration of an elastic beam attached to a cylinder

Following Turek and Hron [49] and Han and Hu [14] a two dimensional FSI benchmark involving flow-induced vibration of an elastic beam attached to a rigid cylinder is considered here. The geometric and setup parameters of the problem are shown in Fig. 13.

The density ratio of the structure to the fluid is $\rho_{s}/\rho_{f}=10$, the Reynolds number $Re=\rho_{f}U_{0}D/\nu=100$, the dimensionless Young’s modulus $E^{*}=E/\rho_{f}U_{0}^{2}=1.4\times 10^{3}$ and the corresponding Poisson ratio $\lambda=0.4$. The initial particle spacing is $dp=0.05D$ and the artificial speed of sound is set to $c_{0}=20{U}_{0}$.

Note that the transport-velocity formulation is applied herein with the background pressure $P_{0}=\rho_{0}c^{2}$. No-slip boundary condition is imposed on the top and bottom walls, while inflow and outflow conditions are applied at the left and right side of the domain, respectively. The inflow condition is given by the parabolic $x$-component velocity profile $U(t,y)=6{U}_{0}y(H-y)/H^{2}$ with a smoothing start technique as in Refs. [49, 14], where $U_{0}=1ms^{-1}$. The inflow and outflow boundary conditions are implemented by the combination of a standard periodic boundary condition and a buffer inflow region with the size of 20 particle spacing by velocity relaxation. After each acoustic time step, within the buffer region, the particle velocity is relaxed to the target value according to the inflow profile by

[TABLE]

where $\mathbf{v}(\mathbf{r}_{i})$ is the prescribe velocity and $\alpha=0.1$ is the relaxation strength chosen here. For the details of the SPH-discretized structure equations, please refer to Ref. [14].

Fig. 14 shows the flow vorticity field and beam deformation at four different time instants in a typical periodic movement when self-sustained oscillation is reached. It can be observed that these results are in quite good agreements with those from previous simulations [49, 14, 50].

Fig. 15 plots the temporal variation of the $x$-, $y$-displacements and the trajectory of the end point $S$, which is defined in Fig. 13. It can be observed that the beam reaches a periodic self-sustained oscillation as time goes beyond a dimensionless time of 50.

A good agreement with computational results of both [50] and [51] is noted. Note that the trajectory of point $S$ shows a typical Lissajous curve with a frequency ratio of $2:1$ between horizontal and vertical components [50].

Table 2 presents a quantitative comparison of the dimensionless amplitude of the oscillation in $y$-direction and the frequency between present computational results and the existing data in the literature [49, 50, 51, 14].

The dimensionless amplitude demonstrates a good agreement, however, slight discrepancies are noted in the dimensionless oscillation frequency. To evaluate the convergence of the proposed method, we performed a convergence study by decreasing the spatial resolution to $dp=0.0125D$. This is found that the amplitude and the frequency, respectively, converge to $0.86$ and $0.186$, again, good agreement with previous results is noted.

4 Concluding remarks

We have proposed a dual-criteria time-stepping method to optimize computational efficiency of the WCSPH method. Apart from the acoustic time criterion which controls the pressure relaxation process, an extra advection time criterion is introduced to control the recreation of particle-interaction configuration. During the advection time criterion, several steps of pressure relaxation are processed whereas the particle-interaction configuration is frozen. As consequence, significant amount of computational efforts are saved due to the fact that particle-interaction configuration can be reused for several time steps. Another noteworthy feature of the present method is that a large CFL number is permissible for pressure-relaxation time-steps within the WCSPH context implying optimized computational performance. Through computational performance analysis, we have demonstrated that an speedup of $2.80$ is achieved by the present method, in comparison to the traditional counterpart. With preliminary examples, we have shown that the present method is capable of modeling free-surface flows involving violent impact and breaking events, as well as fluid- rigid/elastic structure interaction problems, with good robustness and accuracy. While the present simulations are focusing on benchmark tests, the present method is expected to be employed in more complex scientific and engineering problems.

5 Acknowledgement

The authors gratefully acknowledge the financial support by German Research Fundation (Deutsche Forschungsgemeinschaft) DFG HU1527/10-1 and HU1527/12-1 for the present work.