A Simplified Approach of Open Boundary Conditions for the Smoothed Particle Hydrodynamics Method

Thanh Tien Bui,Yoshihisa Fujita and Susumu Nakata

1Graduate School of Information Science and Engineering,Ritsumeikan University,Shiga,525-8577,Japan

2College of Information Science and Engineering,Ritsumeikan University,Shiga,525-8577,Japan

ABSTRACT In this paper,we propose a simplified approach of open boundary conditions for particle-based fluid simulations using the weakly compressible smoothed-particle hydrodynamics (SPH)method.In this scheme, the values of the inflow/outflow particles are calculated as fluid particles or imposed desired values to ensure the appropriate evolution of the flow field instead of using a renormalization process involving the fluid particles.We concentrate on handling the generation of new inflow particles using several simple approaches that contribute to the flow field stability.The advantages of the δ+.-SPH scheme,specifically the particle shifting technique,were successfully applied to correct the position, velocity, and pressure terms of the particles.Therefore, unexpected errors were removed and tensile instabilities of the particles were prevented.The proposed technique is validated for several benchmark test cases,and the tests show that the results match the reference solutions well.A viscous open-channel flow is used to demonstrate the stability of the flow field during the computational time.Based on this stability,we compress the computational domain to a lower resolution in a second test case while preserving the accuracy of the simulation.Flow over a backward-facing step is used to highlight the challenges of inflow boundary conditions with prescribed or non-prescribed values.The developed technique is well suited to the wall boundaries and the evolution of the flow field.The results demonstrate the robustness and versatility of the proposed technique for a variety of simulations.

KEYWORDS Fluid simulation; smoothed particle hydrodynamics; open boundaries

1 Introduction

Open boundary conditions pose an interesting challenge that has attracted a number of investigations by different authors.The simplest and fastest technique, periodic boundary conditions [1,2], is often used in which particles are recycled with the particles passing through the outlet and being re-inserted at the inlet.This technique is limited and is inappropriate for simulations with different inflow/outflow cross-sectional lengths.In addition, violations resulting from the outlet velocity field after re-insertion can cause an unstable state in long runtime simulations.Lastiwka et al.[3] proposed permeable boundary conditions; however, this is difficult to implement with a free surface.A different approach relies on a generalization of the unified semi-analytical boundary condition method introduced by Ferrand et al.[4] and Leroy et al.[5].This technique has been used to impose unsteady open boundaries in incompressible SPH and weakly compressible SPH models.Tafuni et al.[6] used mirrored particles at the fluid zone and a higherorder interpolation process to update the values of the inflow/outflow particles.However, these approaches have a high computational cost and are complicated.Federico et al.[7] presented an implementation of open boundary conditions based on inflow/outflow buffer layers with the inflow particle information assigned at the beginning of the simulation.Alvarado-Rodrıguez et al.[8]utilized a reservoir buffer to prevent a reflection of the velocity field at the outlet.Vacondio et al.[9] introduced open boundary conditions using Riemann invariants.Tarksalooyeh et al.[10]constructed a pre-inlet domain to feed the inflow region.Molteni et al.[11] and Altomare et al.[12] proposed a sponge layer to avoid spurious reflections from waves.Liu et al.[13] used a Taylor series expansion to extrapolate the values of the inlet and outlet particles.Wang et al.[14]and Negi et al.[15] selected a characteristic-based method to apply as a non-reflecting boundary condition in a weakly compressible SPH model.Most previous studies have usually opted to update the inflow/outflow particle information via an extrapolation from the fluid particles.This results in troublesome procedures and is time consuming.

The rapid development of variants of the SPH model, specifically theδ+-SPH scheme [16],has provided reliable results, stability, and high accuracy for a variety of numerical simulations.In theδ+-SPH scheme, density variations were preserved lower than 1% while tensile instability is prevented by applying the particle shifting technique (PST).These advantages of this scheme are used to develop a simplified approach of open boundary conditions in this paper.

To simplify the treatment of the open boundary condition, in the present study, we propose a fast and simple approach based on theδ+-SPH scheme.Instead of implementing a renormalization process using the fluid particles to update the inflow/outflow particle information, we simply treat these particles as the same fluid particles or impose prescribed values to ensure the appropriate evolution of the flow field.This scheme can handle the generation of new inflow particles using several simple approaches combined with the PST.Accordingly, unexpected errors are avoided and the accuracy of the simulations is enhanced.We implemented several simulations to demonstrate the robustness and effectiveness of the proposed technique.

2 The δ+–SPH Model

In the SPH method, the governing equations for weakly compressible SPH in their Lagrangian form are:

whereρ,u,p, µ, andFare the density, velocity, pressure, dynamic coefficient of viscosity, and body force, respectively.

Recently, theδ+-SPH scheme proposed by Sun et al.[16] has emerged as the most effective variant of the SPH model.This scheme is an improvement to theδ-SPH scheme [17,18].The main modification is the use of the PST to remove tensile instabilities, which cause the generation of numerical cavitation.Therefore, this scheme ensures the accuracy and stability of the simulation.Theδ+-SPH model is written as:

whereρi,pi, anduiare the density, ressure, and velocity associated with thei-th particle, respectively,Fgis the body force,Viis thei-th particle volume defined asmi/ρi;Wij.represents the Quintic kernel function andjdenotes the particles inside the support domain of particlei;the gradient ∇iWijindicates the differentiation with respect to the potion ofi-th particle;his the smoothing length;ψijandπijare diffusive and viscous terms;μis the dynamic viscosity;ρ0is a reference density; andc0is the speed of sound.c0.is usually chosen such that:

whereUmaxandpmaxare the maximxpected velocity and pressure, respectively.According to the weakly compressible approach, the density variations must be less than 1%.In addition, the Mach number of the flow should be 0.1 or less.The diffusive and viscous terms are calculated as in[17,18].

We adopted the leapfrog method for the integration time.For stability, the time stepΔt fulfilled the following criteria:

where ‖ai‖ is the particle acceleration and the Courant–Friedrichs–Lewy (CFL)constant was set to CFL=1.5 for all test cases.

The PST was derived by Lind et al.[19] and is indispensable to SPH simulations.In theδ+-SPH scheme, the PST is modified to be well suited for weakly compressible SPH, as shown in Eq.(6).

here, the symbolsrianduiindicate the position and the velocity of thei-th particle, respectively,whiler∗indicates the corrected position after alying particle shiftingδri.Mais the Mach number of the flow.

After the time integration, the pticle position is shifted slightly by an amountδri.to implement the particle shifting algorithm.The corrected position is calculated and obtained at every time step.The pressure and velocity fields are corrected using the Taylor series approximation(see [19]):

whereφis a general variable.To avoid shifting errors, a limit is imposed on the particle shifting distanceδri:

3 A Simplified Approach of Open Boundary Conditions

Open boundary conditions pose an interesting challenge in numerical simulations.The inflow and outflow zones can be defined as being upstream and downstream, respectively, of the fluid domain [7].An upstream instability can seriously influence the development and accuracy of a simulation.The improvement provided by theδ+-SPH scheme allowed us to develop several simple approaches and obtain reliable results [16].

In general, particle-based simulion with open boundary conditions requires four sets of particles that correspond to fluid, wall, inflow and outflow particles as shown in Fig.1.In the methods by Tafuni et al.[6] and Negi et al.[15], the particle motions at inflow/outflow zones are computed using quantities like pressure and velocity of the particles.The particle quantities should be estimated using some kind of extrapolation techniques which requires complicated process like normalized weighted sum of spatial distribution [15] or using ghost nodes for inlet/outlet particles combining with a higher-order interpolation scheme [6].In addition, such extrapolation is required for both of inflow/outflow cases.

Figure 1: Sketch of the simplified approach of open boundary conditions, where kh is the kernel support

In order to avoid the extrapolation process and perform the computation with a simple procedure, we modify the algorithm so that the SPH can be applied not only at the fluid zone but also at inflow/outflow zones.Similar to other inflow/outflow algorithms, we assume that the inflow zone is placed in front of the fluid region so that the attached zone covers a region as wide as the radius of the kernel supportkh, wherehis the smoothing length andkis a parameter varying with the selected kernel function.In our algorithm, the particle motion at the inflow zone is determined as follows:

• The inflow particles move according to their velocity until they cross the border between the inflow and fluid zones.

• The inflow particles that cross the border turn into fluid particles and behave as fluid particles.

• At the same time, the particles that cross the border produce new inflow particles periodically at the front end of the inflow zone and we copy all the information of the particle to the new inflow particle.

Assumingx- andy-axes that are perpendicular and parallel to the border, respectively,the positions of the new inow particles,(xnew,ynew), are determined by(xnew,ynew)=(xinflow−kh,yinflow).In this approach, all the quantities of inflow particles can be determined in the same way as the particles at fluid zone except for the generated particles at the end of the inflow zone.

In terms of the particles at the outflow zone, we apply the SPH algorithm in the same way as the fluid zone and we simply eliminate particles that go out of the outflow zone as shown in Fig.1.Although the total number of particles is not constant because the generation of new inflow particles is not synchronized with the elimination of outflow particles, managing generated/eliminated particles is still simple because the total number of particles does not change so much as far as the domain is fixed.Furthermore, this approach is free from extrapolation process because the quantities of outflow particles is obtained by the SPH algorithm which contributes to simple and consistent computation of fluid particles with open boundary conditions.

In some cases of the extrapolation-based open boundary conditions, prescribed particle quantities given a priori are applied to some of the inflow particles during the entire simulation [6,15].The prescribed values can also be applied in our method with a small modification.In the case where prescribed values are applied to some areas of inflow zone, we simply use the prescribed values instead of copying the quantities of the particles that cross the inflow threshold or applying the SPH algorithm.

4 Numerical Tests

To demonstrate the effectiveness and applicability of the proposed technique, in this section,we perform several numerical test cases.Specifically, we simulate:

– Viscous open-channel flow;

– Flow past a circular cylinder at Re = 200; and

– Flow over a backward-facing step with prescribed and non-prescribed inflow boundary conditions.

4.1 Viscous Open-Channel Flow

A test case of viscous open-channel flow in the laminar regime was conducted to illustrate the stability of the flow over the entire computational domain.The objective of the simulation was to authenticate the steadiness of the velocity field after a sufficiently long period of time and to compare it to the analytical solution shown in Eq.(9).

The initial setup is shown in Fig.2.The fluid is initialized with a fluid domain length ofL=2H, a bottom slope of s0=0.001, a constant density ofρ=1000 kg/m3, and a gravitational acceleration of g=9.81 m/s2, withHbeing the surface depth,μbeing the dynamic viscosity,xandybeing the horizontal and vertical coordinates whose origin is located at the channel bottom.The fluid flow moves with a distribution of velocityu(x,y)for a 2D channel flow, as in Eq.(9), with a null velocity in they-direction.

Figure 2: Computational domain of the open-channel flow

The Reynolds number is calculated such that

where the velocityUis the average velocity profile:

The sound speed was selected to be equal to 10u(y=H), and the dynamic viscosityμcorresponded to the Reynolds number.

The initial boundary condition was imposed as follows:

wherefandidenotes the fluid and inflow particles.The velocity and pressure in the inflow zone were enforced to be the desired values during the entire simulation.At the channel bottom, we enforced a no-slip condition [20,21] using several layers of ghost particles.The simulation was performed for a sufficiently long time to check the stability of the fluid flow.The particles were initially distributed with a resolution of 4Δx=4H/125 (see [7])with approximately 2000 particles with Reynolds numbers ofRe=10 and 100.

Figs.3a and 3b illustrate the particle distribution and velocity field att(g/H)1/2=100 forRe=10 and 100, respectively.The obtained velocity fields are smooth and provide good results without numerical noise or instabilities.The velocity vector field in Fig.3c shows that the fluid develops with parallel layers over the entire comtational domain.

Figure 3: Particle distribution and velocity fields at t(g/H)1/2=100 for (a) Re=10 and (b) Re=100, and (c)velocity vector field for Re=100

To verify the stability, three differentx-positions,x=0 (inflow threshold),x=H(middle of the fluid domain), andx=2H(outflow threshold), were considered to compare the analytical solution in Eq.(9)to the numerical results obtained via the SPH.Fig.4 shows the good agreement between the analytical and numerical results.The obtained results give approximately the same analytical solution at all three differentx-positions.This demotrates the steadiness of this technique throughout the entire computational domain.

Figure 4: Comparisons between the analytical solution and the numerical results at t(g/H)1/2=100 for Re=10 and 100 at x=0, x=H, and x=2H

To check the convergence of the velocity field between the proposed technique and the analytical solution, the mean square error percent (MSEP)was calculated using Eq.(13)atx=H(in the middle of the fluid domain)with a resolution of 4Δx:

whereuaandunare the analytical and numerical velocities, respectively, andNis the number of velocity values.

Fig.5 illustrates a comparison of the MSEP values forRe=10.In the figure, the MSEP value of the proposed technique is low with a maximum error of approximately 0.1% at a resolution of 4Δx, while that of Federico et al.[7] has a maximum error of 3.2%.The MSEP values are significantly lower than with those of Federico et al.[7], even for a resolution ofΔx.Fig.6 depicts the MSEP value of the velocity field of the proposed technique atRe=100, howing a peak error of approximately 0.12%.

Figure 5: The mean square error percent (MSEP)for Re=10 showing the proposed technique and the inflow/outflow modeled by Federico et al.[7]

Figure 6: MSEP for Re=100 using the proposed technique

4.2 Flow Past a Circular Cylinder at Re=200

The second test case consistedf flow past a circular cylinder, where the circular cylinder is represented by an implicit function:

whereris the radius of the circular cylinder.The implicit function of a circular cylinder is defined as the zero-level set of the functionf(x,y)=0 and is placed inside a rectangular channel.The boundary treatment technique in [22] was applied to generate wall particles.The boundary conditions were implemented as follows:

(1)No-slip boundary conditions were enforced on the body surfaces, and the information concerning these solid particles was extrapolated from the fluid particles following the method in [21].

(2)For the upper and lower walls, symmetry boundary conditions for the velocity were imposed as described in [23], i.e.,vy=0 and∂vx/∂y=0.

Fig.7 depicts the initial setup of the simulation.A smaller domain and lower resolution compared to another study [6] were considered to prove the stability and effectiveness of the proposed technique.Following Lastiwka et al.[3], the computational domain was [0,15D]×[0,8D]with a cylinder with a diameter of D=0.1 m located at the position [5D,4D].A constant density ofρ∞=1000 kg/m3, a constant x-velocity ofU=1 m/s, and a null y-velocity were initialized for the fluid particles.The fluid phase was discretized with the initial particles placed at a resolution of D/Δx=40.The quintic kernel function was used with the smoothing length set toh=1.5Δx.The Reynolds number of the flow wasRe=UD/ν, whereνis the kinematic viscosity in m2/s.

Figure 7: Computational domain for 2D flow past a circular cylinder

To prevent an impulsive start of the velocity field, a constant accelerationa0=U2/Dwas used fort

The drag and lift force coefficients on the body are calculated such that:

whereFxandFyrepresent the total horizontal and vertical forces on the body, respectively,following the method in [25].

Validation for this case was performed by comparing the time histories of the drag and lift coefficients using the technique proposed in Tafuni et al.[6], as shown in Fig.8.Good agreement is obtained between the two solvers.The magnitude of the drag force increases to a peak value and then decreases to a steady state during the time intervaltU/D∈[0,40] atRe=200.The drag and lift coefficients converge to a steady values ofCD=1.46 andCL=±0.73, respectively, forRe=200.These results are similar to those obtained by other authors, such as Marrone et al.[25],Vacondio et al.[26], and Negi et al.[15].

Figure 8: Time history of the drag and lift coefficients for flow past a circular cylinder at Re=200

Fig.9 illustrates the velocity and pressure contours atRe=200.Two dimensionless values, the velocity magnitudeU∗=u(x/D,y/D)U−1andP∗=2P(x/D,y/D)ρ−1U−2, are used.The velocity and pressure fields are smooth throughout the entire computational domain.The tensile instability that causes the generation of numerical cavitation was eliminated by the PST.Fig.9 depicts an unsteady von Karman street with an oscillatory wake that appears behind the cylinder atRe=200.

Figure 9: 2D flow past a cylinder: the velocity and pressure fields at Re=200

Fig.10 shows the velocity vector field atRe=200.It can be seen that the velocity vector field at [0,3.5D] on the horizontal axis developed with parallel layers and exhibits stability without numerical noise.This result was obtained as aesult of the upstream stability, where the proposed scheme was applied.

Figure 10: 2D flow past a cylinder: the velocity vector field at Re=200

4.3 Flow over a Backward-Facing Step

4.3.1 2D Backward-Facing Step with Prescribed Inflow Boundary Conditions

A 2D backward-facing step problem was simulated atRe=389 and compared to the experimental results of Negi et al.[15] and Armaly et al.[27].Following Negi et al.[15], the inflow particles were set to a constant velocity for the entire computational time.However, no-slip boundary conditions were enforced at the walls, causing non-physical pressure fluctuations.To address this issue, a small part of the initial wall next to the inlet zone can include an imposed slip, as in [15], or a parabolic profile of the velocity can be used at the inlet zone, as in [8,28].These approaches can cause several inconveniences during test cases.Taking advantage of the proposed technique, we only imposed a prescribed velocity value at the initial time step and then the velocity automatically adjusted to be well suited to the near-by boundary conditions.The velocity in the inflow zone was less than the maximum velocity and ensured a Mach number of the flow ofMa<0.1.

The step height was set toH=4.9 mm with an inlet width ofH1=5.2 mm.The density wasρ=1.0 kg/m3, the dynamic viscosity wasμ=0.017 kg/ms, and the length of the domain wasL=130.6 mm.The initial inlet velocity was 0.67 m/s with a maximum of 1 m/s.We compared the velocity profiles at four different marked positions P1–P4 ofx/H= 2.55, 3.57, 4.80, and 7.14(where x is the distance downstream of the step).A schematic of the simulation model is shown in Fig.11.

Figure 11: Sketch of the domain used for the backward-facing step simulations

Fig.12 illustrates the axial velocity of the simulation atRe= 389.The maximum axial velocity is 1.0 m/s, and the minimum value is −0.11 m/s.The PST helps remove particle clumping near the step.Good performance for the axial velocity is shown, where the velocity at the inlet automatically adjusts corresponding to the wall boundary with the velocity profile shown in Fig.13a.The velocity tends to low near the wall and to reach maximum values at the center of the domain.This is likely equivalent to using a special treatment at the inlet, as in [15].

Figure 12: Axial velocity for Re=389

Figure 13: (a)Velocity profile at the inlet, (b)reattachment length for Re= 389, (c)Time history of number of particles in the entire computational process

The reattachment length was determined to bex/H=7.84, as shown in Fig.13b.This result is close to the result obtained by Negi et al.[15] of 7.9 and the experiment value [27] of 7.94.

Fig.13c shows the variation of fluid particles in the entire computational process.The simulation starts with 44550 particles and reaches to steady state after around 120000 steps at which the number of the particles is 43599 particles.The percentage of mass error is 2.13%.

The velocities at the four different marked positions, P1–P4, were considered over the channel height and compared to the profiles of the reference results from [15,27], see Fig.14.The results obtained are in good agreement with the reference solution at each marked position.

Figure 14: Comparison of velocities at four different locations for Re=389

4.3.2 2D Backward-Facing Step with Non-Prescribed Inflow Boundary Conditions

In this case, a 2D backward-facing step problem was simulated atRe=100 and compared to the experimental results of Adami et al.[21] and Issa et al.[29].In these studies, another step was added at the end of the channel with the same length to apply periodic conditions.Using the proposed technique, we can easily simulate this case without a second step.

In this simulation, the density wasρ=1.0 kg/m3, the viscosity wasν=1.456×10−2, and the length of the domain wasL=122.4 m.We applied a constant body force ofF=3.0×10−4m/s and a sound speed ofc=2.1.The step height was set toH=4.9 m with an inlet width ofH1=5.2 m.No-slip boundary conditions were enforced on the walls.A schematic of the simulation model is shown in Fig.15.We compared the velocity profiles at four different marked positions,P1–P4, as shown in Fig.15.

Figure 15: Sketch of the domain used for the backward-facing step simulations

The fluid flow was driven by a constant body force F.The fluid particles propagated and evolved during the simulation until reaching a corresponding stable state atRe=100.This led to unknown values in the inflow region.The non-prescribed inflow test case demonstrates the effectiveness and applicability of the proposed technique.

Fig.16 illustrates the axial velocy of the simulation atRe=100.It can be seen that the flow field develops smoothly over the entire mputational domain.The fluid particles are colored with the maximum axial velocity being 0.21 m/s and the minimum value being −0.015 m/s.

Figure 16: Axial velocity for Re=100

Assessing the quality of simulation, the axial velocities at four different marked positions, P1–P4, were considered over the channel height and compared to the reference profiles from [21,29],as shown in Fig.17.The results obtained at each marked position are in fairly good agreement with the reference solution.

Figure 17: Comparison of the axial velocities at four different marked positions, P1–P4, for Re=100

Fig.18 depicts the variation of fluid particles in the computational time.The number of fluid particles is 26160 at the beginning and it reaches to 26281 particles atRe=100.The percentage of mass error is 0.46%.

5 Conclusions

In this paper, we developed a simplified approach of open boundary conditions.Taking advantage of theδ+-SPH scheme, the proposed technique is simple to implement but is effective and robust in test cases.The approach is classified into sets of fluid, wall, inflow, and outflow particles, as in most inflow/outflow algorithms.The inflow/outflow information is treated as fluid particles to ensure appropriate evolution to the flow field instead of using a renormalization process referencing the fluid particles.The values of new inflow particles are copied from the shifted particles, and the differences are corrected using the PST.

The numerical tests demonstrate that the proposed technique obtains good results with a high agreement with the reference solutions.The first test case demonstrated the stability of the flow field over a sufficiently long time with the MSEP value of the velocity field being approximately 0.1%.Given this stability, we compressed the computational domain to a lower resolution in a second test case to demonstrate the high accuracy of the simulation.The third test case examined the flexibility of the inflow boundary conditions to prescribed or non-prescribed values.This case developed results well suited to the wall boundary and the evolution of the flow field.The performance of the results demonstrated the versatility of the proposed technique.Overall, the proposed technique has addressed the information on inflow/outflow particles without extrapolation from the fluid zone.

One of the limitations is that particles in the fluid zone can move across the borders of inflow zones in case of flow reversion.This issue is pointed out also in [6] and a decent solution is given in the paper.We expect to improve this issue in future work.

Acknowledgement:This research was supported by the Ministry of Education, Culture, Sports,Science and Technology (MEXT), Japan.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.