- Email: [email protected]
On information coupling in hybrid ISPH framework for fluidized granular systems
Advances in Water Resources 131 (2019) 103366
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Short communication
On information coupling in hybrid ISPH framework for fluidized granular systems Gourabananda Pahar a,∗, Anirban Dhar b a b
Assistant Professor, Department of Civil Engineering, Indian Institute of Technology Kanpur, Kanpur UP 208016, India Associate Professor, Department of Civil Engineering, Indian Institute of Technology Kharagpur, Kharagpur WB 721302, India
a r t i c l e
i n f o
a b s t r a c t
Keywords: ISPH Sediment transport Information coupling
An improved information transfer mechanism is developed to couple fluid and granular modules for Incompressible Smoothed Particle Hydrodynamics modelling of sediment transport. The interaction force pair necessitates effective projection of information among parallel continuum based modules. The particle position in either module is susceptible to considerable variation due to their Lagrangian nature. An attempt has been made to minimize the effective thickness of the diffused interface between the fluid and granular particles. The proposed model reduces artificial drag force resulting from diffused computation of effective porosity. The revised framework outperforms in coarse resolution compared to its existing counterpart.
1. Introduction
the fraction of the volume associated with each fluid cell. Traditionally CFD-DEM assumes the full volume of the solid particle stays within the corresponding grid. Thus, all solid particles present within a fluid cell have equal contribution through their representative volume (irrespective of their positions). Peng et al. (2016) have introduced a Particle Meshing Method to take care of solid particle volume present within an Eulerian grid. In a Lagrangian method like SPH, the fluid cell volume can be mapped to the idea of influence domain. The contribution of particles has been generally taken through an uncorrected kernel value employed on them (Sun et al., 2013). The contribution of the particle decreases monotonically from the concerned location leading to an efficient representation. Typically, CFD-DEM is an ensemble of coupled Eulerian (CFD) and discrete (DEM) methods. The effective porosity and other coupling parameters are directly (or crudely) interpolated at the center of the fluid cell (Sun and Xiao, 2016; Zhao et al., 2017). Thus the resolution of fluid cells inevitably becomes coarser compared to the discrete granular module. In coupled SPH-DEM (Robinson et al., 2014; Markauskas et al., 2017; Tan and Chen, 2017) systems, particles in both parallel modules are moving separately. Thus the information needs to be suitably interchanged between modules. The drag force computation for any particle requires relative velocity at a common location, which is computed generally through a Shepard corrected kernel (Robinson et al., 2014). The pressure gradient and viscous term acting on the solid particle can be evaluated directly for a coupled CFD-DEM method from the corresponding fluid cell. Coupled SPH-DEM (Markauskas et al., 2017; Tan and Chen, 2017) generally avoids the pressure coupling by directly incorporating buoyancy force into account. It also removes the
Lagrangian particle methods are generally well-equipped to capture the interaction between rapid flow and fluidized granular bed (e.g., landslide-generated waves, sediment transport). Coupled flow models utilize various algorithms for interchanging information between submodules (.i.e. fluid and granular media). However, sharp discontinuity is often flattened out by these interpolation algorithms. This work aims to reduce the adverse impact of information loss by improving interpolation algorithms between different components of the coupled flow framework. The continuum conceptualization of effective porosity does not remain valid beyond a certain minimal representative elemental volume (REV). The volume-averaged flow equations lose their validity beyond this threshold length-scale corresponding to minimum REV. The distribution of effective porosity across interface should not depend on the resolution adopted for a particular simulation. Several researchers have formulated various elaborate algorithms related to Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) to achieve an actual distribution of effective porosity. Wu et al. (2009) have developed an efficient analytical method based on particle volume shared by various unstructured grids for coupled discrete particle methods. Generally, CFD-DEM directly calculates voidage function through the number of particles present in a particular fluid cell (Zhao and Shan, 2013; Shan and Zhao, 2014). Recent studies have suggested creating a timedependent mesh/ analytical solution (Boyce et al., 2014) based on current co-ordinates of granular particles at every time step to calculate
∗
Corresponding author. E-mail addresses: [email protected] (G. Pahar), [email protected] (A. Dhar).
https://doi.org/10.1016/j.advwatres.2019.06.011 Received 4 July 2018; Received in revised form 1 May 2019; Accepted 24 June 2019 Available online 26 June 2019 0309-1708/© 2019 Elsevier Ltd. All rights reserved.
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
likely to satisfy as particles tend to travel separately in their respective domains. Unlike the existing SPH-DEM formulations, the present framework utilizes initial identical (and coarser) resolution for the granular module. The improved technique utilizes a first-order consistent method to take care of rapid yet effective interpolation of information between parallel modules. Effect of artificial drag forces resulting in the previous model has been thoroughly minimized by the proposed variant. The modified framework is compared with existing experimental data along with numerical result simulated by the erstwhile formulation. 2. Information transfer between modules Anderson and Jackson (1967) had developed the averaged equations for fluid flow over deformable granular bed. Considering continuum approximation for the granular phase, the governing equations of continuity and momentum for granular and fluid phases can be written as (Pahar and Dhar, 2018), Fig. 1. Coupled modular structure.
possibility of enforcing the erroneous pressure gradient resulting from spurious numerical fluctuations in weakly compressible SPH. Generally, locally averaged Navier-Stokes equation requires a comparatively coarser resolution (compared to the granular module) of fluid phase to obtain sufficiently smooth porosity field (Sun et al., 2013; Tan and Chen, 2017). There are also instances of coupled SPH-DEM frameworks (Komoroczi et al., 2013; Ren et al., 2014; Canelas et al., 2016), where DEM resolution is coarser compared to SPH. However, they are not suitable for simulation of sediment-transport mechanism. A hybrid continuum based ISPH-ISPH framework has been developed in a previous study (Pahar and Dhar, 2017) for incompressible modelling of sediment transport. The framework shares similarities with existing SPH-DEM models though the granular media is modelled using an incompressible continuum conceptualization (Fig. 1) through rheological approximation. The modules interact with each other through an interaction force pair (Kafui et al., 2002). The force pair consists of fluid pressure-gradient, the divergence of viscous stress tensor and drag force. Each of these components would require extensive information from the respective parallel module. The fluid pressure is inherently smooth due to divergence-free nature of truly incompressible SPH (Lee et al., 2008). Thus fluid pressure gradient is incorporated into the interaction force pair without the additional assumption of hydrostatic pressure distribution (buoyancy). Different shades of blue in the continuum based fluid phase in Fig. 1 represent flow through porous media (lighter) and outside it (darker). The interface between two phases should be clearly delineated to distinguish the effect of resistive force vector originating from solid skeletal structures. The modules exist in parallel domains although their coordinate systems coexist. Interaction forces are calculated based on the projection of one module to another one. Prior work utilized a zeroth order smoothing algorithm to transmit information between modules. Though It ensures a smoothed distribution of different properties, the rapid discontinuity is flattened out leading to loss of information. For example, drag force varies non-linearly with corresponding effective porosity. It vanishes for zones having unit effective porosity. A large diffused interface between porous and non-porous zone would exert an extra artificial drag force on the fluid module. Thus it may adversely affect the computation of erosive processes. In the current study, an improved information transmission mechanism between coupled modules has been proposed. The projection of information becomes tricky for these separated yet coupled modules as both sets of particles are highly susceptible to irregular distribution. Gradient/ divergence(s) of properties are directly considered for projection in the higher order interpolation. SPH requires information regarding volume associated with the center particle to calculate the gradient/ divergence at the concerned location. However, these criteria are un-
∇ ⋅ 𝐮𝑠 = 0
(1a)
∇ ⋅ (𝜂𝐮𝑓 ) = 0
(1b)
𝐅𝑓 𝑠 𝐷𝐮𝑠 1 1 =− ∇𝑃𝑠 + 𝐠 + ∇ ⋅ 𝛕𝐬 + 𝐷𝑡 𝜙𝜌𝑠 𝜙𝜌𝑠 𝜙𝜌𝑠
(1c)
𝐷(𝜂𝐮𝑓 ) 𝐷𝑡
=−
𝐅 𝜂 1 ∇𝑃𝑓 + 𝜂𝐠 + 𝜈∇2 (𝜂𝐮𝑓 ) − 𝐷 + ∇ ⋅ 𝛕̄ 𝜌𝑓 𝜌𝑓 𝜌𝑓
(1d)
where 𝐮𝑠 = granular velocity, 𝐮𝑓 = fluid velocity, 𝜙 = solid fraction, 𝑃𝑠 = granular pressure, 𝐠 = acceleration due to gravity, 𝜌𝑠 = granular density, 𝐅𝑓 𝑠 = interaction force pair per unit volume, 𝜂 = effective porosity, 𝑃𝑓 = ̄ fluid pressure, 𝐅𝐷 = drag force per unit volume, 𝜈 = fluid viscosity, 𝜌𝛕 = 𝑓
sub-particle-scale tensor for modelling small-scale flow information. The interaction force can be represented as Kafui et al. (2002), 𝐅𝑓 𝑠 = −(1 − 𝜂)∇𝑃𝑓 + (1 − 𝜂)∇ ⋅ 𝛕 + 𝐅𝐷
(2)
Drag force (FD ) is determined through the following expression following Di Felice (1994). 𝐅𝐷 =
1 𝑛 𝐶 𝜌 𝜋𝑑 2 𝜂 (2−𝜉) (𝐮𝑓 − 𝐮𝑠 )|𝐮𝑓 − 𝐮𝑠 | 8 𝑝 𝑑𝑖 𝑓
(3)
where 𝑛𝑝 = number of granular particles present in a unit volume, 𝜉 = corrective porosity function. The drag force component is a non-linear function of effective porosity. In a continuum approach, np can be calculated as the volume of solid (1 − 𝜂) divided by the volume of a single particle ( 𝜋6 𝑑 3 ). Generally, bed-roughness is treated through drag-based formulation (Kazemi et al., 2017) near the bed-zone. The bed is considered as a rigid boundary for these models. The proposed model does not consider granular media as a solid boundary. Rather fluid outside and inside porous media is solved through a volume-averaged filtered Navier-Stokes equation. The momentum equation contains drag force which incorporates the information regarding particle diameter and relative velocity between fluid and granular particles. A diffused interface will ensure a smudged distribution of np , which in turn introduces artificial drag force even when the granular particles are not immediate neighbours of fluid particles. The information transfer mechanism between modules is described in detail in the following sections. 2.1. Information transfer: granular to fluid module The terms containing fluid pressure gradient and divergence of viscous stress are implicitly embedded into the fluid module. Thus, the informations required from the granular module for calculating drag force are the effective porosity and granular velocity at each fluid particle. According to Taylor Series, any quantity 𝜓 can be written as, | 𝜓𝑏 = 𝜓𝑎 + 𝐫𝑎𝑏 ⋅ ∇𝜓 | + ⋯ (4) |𝑎
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
where 𝐫𝑎𝑏 = relative position vector between particles a and b, 𝐫𝑎 = position vector of particle a. In discrete form the aforementioned equation can be written as, 𝜓𝑏 = 𝜓𝑎 + 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑ (𝜓𝑐 − 𝜓𝑎 )𝑣𝑐 ∇𝑤𝑎𝑐
(5)
𝑐=1
𝑐=1
Eq. (6) is valid for only a pair of particles inside the influence region. Similar sets of equations can be produced by sweeping over all the particles present inside support domain. Applying a kernel throughout the particles inside the support domain, the equation can be modified as, 𝑁𝑡 ∑
[(
𝑏=1 𝑁𝑡
=
∑ 𝑏=1
1 − 𝐫𝑎𝑏 ⋅ [(
| 𝜙𝑏 = 𝜙𝑎 + 𝐫𝑎𝑏 ⋅ ∇𝜙| |𝑎 = 𝜙𝑎 + 𝐫𝑎𝑏 ⋅
𝑐=1
𝑣𝑐 = volume associated with particle c present in the support domain of particle a, 𝑤𝑎𝑐 = kernel value at any particle c with respect to center particle a. The total number of particles present in a compact support domain (of a single module) is denoted by Nt . The center particle position is represented by a. However, the information regarding center particle position and volume associated with it remains unknown. All the particles present in the support domain except the concerned center particle remain in an identical module. The kernel gradient possesses zero value at the center location. Thus, Eq. (5) can be written as, ( ) ( ) 𝑁𝑡 𝑁𝑡 ∑ ∑ 𝜓𝑎 1 − 𝐫𝑎𝑏 ⋅ 𝑣𝑐 ∇𝑤𝑎𝑐 = 𝜓𝑏 − 𝐫𝑎𝑏 ⋅ 𝜓𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 (6)
𝜓𝑎
the domain, the solid fraction at central location can be written as,
𝑁𝑡 ∑ 𝑐=1
𝜓𝑏 − 𝐫𝑎𝑏 ⋅
)
]
𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏 𝑁𝑡
∑ 𝑐=1
)
(7)
2.1.1. Effective porosity calculation In CFD-DEM, effective porosity of a fluid cell is generally determined based on the total volume of solid particles present in the respective cell. To obtain the same in the coupled SPH-SPH model, a concept of ghost particles is introduced (Fig. 2). Fig. 2a depicts the representative particle positions for both fluid (blue) and granular (yellow) modules. The presence of granular particles in the support domain has been shown in Fig. 2b. A discontinuity of solid fraction exists exactly at the interface. The completeness of the corrupted support domain is maintained using a few virtual (grey) particles. The value of solid fraction associated with these particles is zero. The gradient (and higher order derivatives) of solid fraction is also considered to be zero. The location and associated volume of these ghost particles remain unknown. The representative volume associated with virtual granular particles need not be identical to that of actual granular particle. Outside the porous media, For the real particles present inside
Fig. 2. Effective porosity calculation.
𝑐=1 𝑁𝑡
∑ 𝑐=1
(𝜙𝑐 − 𝜙𝑎 )𝑣𝑐 ∇𝑤𝑎𝑐
(8)
𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐
Nt represents the total number of particles present in the compact support domain. Thus Nt becomes equal to 𝑁𝑎𝑐𝑡 + 𝑁𝑣𝑖𝑟 , where Nact and Nvir are respective numbers of actual and virtual particles inside the support domain. −𝜙𝑎 can be discarded as the influence domain has compact support (zeroth order consistency). For the virtual particles assumed inside the domain, the solid fraction at central location can be written as, | 𝜙𝑎 = 𝜙𝑏 + 𝐫𝑎𝑏 ⋅ ∇𝜙| |𝑏
(9)
Eq. (9) utilizes the information regarding the solid fraction and gradient of solid fraction (which are zero as discontinuity prevails outside granular domain) at the virtual particles. Accordingly, these particles have zero contribution at central location. The ghost particles also en∑𝑁 sures 𝑏 𝑡 𝑤𝑎𝑏 𝑣𝑏 to become unity for compact support domain. Thus solid fraction at any location can be written as, 𝜙𝑎 =
]
𝜓𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏
= 𝜙𝑎 + 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑
𝑁𝑎𝑐𝑡
∑
[(
𝑏=1
𝜙𝑏 − 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑ 𝑐=1
)
]
𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏
) ] ∑ [( | 𝜙𝑏 + 𝐫𝑎𝑏 ⋅ ∇𝜙| 𝑣𝑏 𝑤𝑎𝑏 |𝑏 𝑏=1 [( ) ] 𝑁𝑎𝑐𝑡 𝑁𝑎𝑐𝑡 ∑ ∑ 𝜙𝑏 − 𝐫𝑎𝑏 ⋅ 𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏 = 𝑁𝑣𝑖𝑟
+
𝑏=1
(10)
𝑐=1
For the virtual granular particles, solid fraction and gradient of solid fraction are zero. It may appear that zero solid-fraction and its gradient for virtual granular particles is rather conflicting. However, the calculation of gradient depends upon choice of smoothing length and representative volume associated with virtual particles. Considering very small smoothing length along with numerous virtual granular particles represent physically consistent result. Their number, location and associated
Fig. 3. Effective porosity distribution.
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
Fig. 4. Particle distribution.
Fig. 7. Information transfer: from fluid to granular media.
Fig. 5. Effective porosity distribution: (a) Direct averaging, (b) Direct averaging with edge trimming, (c) Zeroth order, (d) First Order, (e) actual.
Fig. 6. Error in calculation of effective porosity distribution: (a) Direct averaging, (b) Direct averaging with edge trimming, (c) Zeroth order, (d) First Order.
volumes do not feature in the equation. Thus, effective porosity at any point can be written as, 𝜂𝑎 = 1 − 𝜙𝑎
(11)
2.1.2. Velocity interpolation Granular velocity also needs to be interpolated on the fluid particle position for relative velocity calculation. Including the virtual particles in the simulation would require knowledge of their respective velocity, which is likely to be non-zero (unlike the previous case). Thus, the virtual particles cannot be directly included. The kernel gradient has to be corrected for corrupted support domain. The interpolated property 𝜓
Fig. 8. Initial conditions.
can be written as (Eq. (7)), ) ] ∑𝑁𝑎𝑐𝑡 [( ∑𝑁𝑎𝑐𝑡 𝜓𝑐 𝑣𝑐 ∇𝑤̃ 𝑎𝑐 𝑣𝑏 𝑤̃ 𝑎𝑏 𝜓𝑏 − 𝐫𝑎𝑏 ⋅ 𝑐=1 𝑏=1 𝜓𝑎 = ∑ [( ) ] ∑𝑁𝑎𝑐𝑡 𝑁𝑎𝑐𝑡 𝑣𝑐 ∇𝑤̃ 𝑎𝑐 𝑣𝑏 𝑤̃ 𝑎𝑏 1 − 𝐫𝑎𝑏 ⋅ 𝑐=1 𝑏=1
(12)
𝑤̃ and ∇𝑤̃ represent Shepard corrected kernel and corrected (not Shepard corrected) kernel gradient (Bonet and Lok, 1999). The kernel gradients do not need to be zeroth order consistent, as the difference
G. Pahar and A. Dhar
(𝜓𝑐 − 𝜓𝑎 ) is incorporated for computation of gradients. Zeroth and firstorder consistent kernel gradients will contain volume and location of the unknown central particle. Utilizing first order consistent kernels will implicitly interpolate the information of central particle. Representative distributions of effective porosity for different methods have been shown in Fig. 3. Fig. 3a depicts actual distribution of solid particles and corresponding effective porosity associated with fluid particles. The porosity can be calculated directly by averaging the void space in a unit volume (Fig. 3b). However, it leads to a diffused interface having width 𝜍, where 𝜍 denotes the diameter of the circular support domain. The partition of volume associated with fluid particles inside the support domain affects the calculation of effective porosity. The inverse distance weighting property of the kernel reduces the diffused nature of the interface (Zeroth order). Including the kernel gradients into the Eq. (7) would lead to a more rapid transition of effective porosity. Ideally, an abrupt change of effective porosity should exist between these zones. However, in a continuum approach, it is difficult to calculate the actual distribution of discontinuous effective porosity. Effective porosity cannot be calculated by utilizing the nearest neighbour approach (with minimal average in-
Advances in Water Resources 131 (2019) 103366
terparticle distance) due to its inherent continuity and differentiability constraints. To validate the proposed algorithm a representative coupled model has been shown. Distributions of fluid and granular particles have been depicted in Fig. 4. A regular particle distribution has been adopted for both modules with identical inter-particle distance. The topmost granular particle layer superimposes with the fluid particle layer. The distributions of effective porosity obtained from direct averaging (with and without volume trimming), kernel and kernel gradients have been shown in Fig. 5a–d. The exact distribution of effective porosity for the fluid particles has been depicted in Fig. 5e. The effective porosity of fluid particle layer coinciding with topmost granular layer is considered to be (1 − 𝜙2 ). Direct averaging considers identical contributions of particles present in a support domain. Each particle is considered to be associated with an identical circular volume (Fig. 5a). The effective porosity is calculated from the ratio of the total volume represented by the solid particles and the volume of the support domain (following traditional CFD-DEM methods). It underestimates the effective porosity for the interior portion of the granular domain. Total volume carried by every
Fig. 9. Effective porosity distribution at, (a) t = 0.25 s, (b) t = 0.50 s, (c) t = 0.75 s, (d) t = 1.00 s, (e) t = 1.25 s, (d) t = 1.50 s .
G. Pahar and A. Dhar
particle is included irrespective of its position. Total volume associated with the particles present in the vicinity of the boundary of the support domain may not remain within the support domain. Edges are trimmed for the second case (Fig. 5b). However, the width of the interface remains just as large as the previous case. The circular volume representation becomes essentially invalid for the current particle distribution. Each particle represents a square volume with the particle placed at the center of the square. A Voronoi diagram should be drawn at each timeinstance around the particles to get accurate volume representation as particle position is susceptible to large temporal variation. A straightforward kernel approximation is adopted in the third case (Fig. 5c). The corruption in the support domain is calculated through completeness of kernel. However, particle layers present above granular profile get heavily influenced. The current method provides suitable distribution of effective porosity leading to sharp interface (Fig. 5d). The error (Δ𝜂) in effective porosity calculation for each of these methods have been shown in Fig. 6, where Δ𝜂 = 𝜂exact − 𝜂calculated . The embedded kernel gradient influences the effective porosity through corresponding particle distance concerning the center particle. The kernel gradient peaks in the vicinity of the center particle (for both compact/ corrupted support domain). In SPH, a continuous kernel is transformed into discrete quantities for weight calculation. Thus a larger support domain would be a better choice regarding unity representation. However, the width of the diffused interface also increases with larger support domain. In the current study, 1.60Δx has been taken as smoothing length for a C2 kernel (Wendland, 1995). Maximum 6 layers of particles will be influenced by the proposed method for identical resolution of either module. The effect of kernel gradient vanishes for a full influence domain (the interior portion of granular media). The error in the calculation of effective porosity for a compact support domain is in order of 0.47%. Adopting higher order kernels (C4, C6) may further decrease the error in computation. However, being too concen-
Advances in Water Resources 131 (2019) 103366
trated around the central particle, they would require a larger support domain effectively increasing the width of the diffused interface. It should be also noted that the mechanism (for relative velocity calculation) is developed on the basis of zeroth order inconsistency ∑ 𝜕 𝑤̃ 𝑎𝑐 ( 𝑁 𝑐=1 𝜕𝑥 ≠ 0). Both Shepard corrected kernel and kernel gradient must include the center point and volume associated with it. In the current coupled modular framework, fluid and granular particles are present in the parallel domain while projecting information to another counterpart. Thus the volume of the projected center particle remains unknown. The volume of the center particle does not feature in the calculation of the proposed method. 2.2. Fluid to granular module Unlike the last case, the full interaction pair needs to be incorporated into the granular module. Thus the pressure gradient and divergence of stress have to be accurately calculated at every granular particle position. The procedure has been shown in the Fig. 7. The transfer process has been highlighted for a few isolated granular particles present inside the high-velocity fluid regime. Generally, the granular domain remains fully submerged leading to a compact support domain (in parallel modules). For erosion resulting from the highvelocity fluid stream, granular particles may become separated from the bed. The fluid property is calculated on the granular particles with zeroth order averaging through Shepard corrected kernel in the previous framework. Ultimately, the gradient/divergence of any property can be computed based on the support domain of the granular media in the final step. This procedure will be able to produce proper result inside the porous media; i.e., where each granular particle has a sufficient number of granular particles present in its support domain. However, for isolated granular particles, this may result in erroneous interpolation of information. The case shown in Fig. 7 contains 3–4 particles present
Fig. 10. Effective porosity distribution at, (a) t = 0.25 s, (b) t = 0.50 s, (c) t = 0.75 s, (d) t = 1.00 s for Scenarios-IIIa (identical resolution) and IIIb (coarse granular resolution).
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
divergence on each fluid particles. However, the procedure is completed within the first two steps of the aforementioned mechanism. 3. Results and discussions
Fig. 11. Initial distribution of effective porosity (𝜂) for Scenario-III.
in the influence zone each granular particle. Sparse irregular distribution of these particles (around the edge of the support domain) will result in a badly-scaled corrective matrix for kernel gradient modification (Bonet and Lok, 1999). Thus the kernel gradient may have to be taken as zero for these secluded particles. However, inside a highvelocity fluid stream, the pressure gradient is bound to be sufficiently high to affect the particle movement. Thus, a different procedure has been proposed in the current study. In the present work, the pressure is not directly interpolated into the granular particles. Instead, the fluid pressure gradient is calculated on each fluid particle. The pressure gradient is directly interpolated on the granular particles. The procedure requires additional computational overhead for searching and calculating pressure gradient/stress
The proposed mechanism is applied to three cases of dam-break induced erosion enclosed in a box. The initial condition is depicted in the Fig. 8. All three cases have sand as the bed material (Spinewine, 2005). The lock-gate is removed swiftly to study the temporal variation of bed. These cases are also simulated using the previous method and compared with the proposed variant. The resolution is deliberately kept coarse (initial interparticle spacing l0 = 0.0125 m) to assess the effectiveness of the model. All the simulation parameters (e.g. support domain radius, free-surface detection criteria, etc) have been kept identical for both variants. The third case is also replicated with separate resolutions for different modules (𝑙0 |𝑠𝑎𝑛𝑑 = 0.02𝑚, 𝑙0 |𝑤𝑎𝑡𝑒𝑟 = 0.0125𝑚) to study the effect of the developed algorithm. The smoothing length is taken as 1.60l0 . Calculation of effective porosity utilizes the smoothing length of the granular domain, whereas the transfer of fluid information to granular module considers smoothing length of the fluid module. In case of identical resolutions, the smoothing length for either of the cases remains unchanged. Simulated fluid profiles with the proposed framework along with experimental data have been shown in Fig. 9. The effective porosity serves as an indicator of the granular profile at different time instances. The experimental data of Spinewine (2005) has been shown by using circles (for water) and triangles (for granular profile). The progression of fluid fronts for the third scenario along with identical/ coarse resolutions (for granular media) has been depicted in Fig. 10 The coupled ISPH framework with the proposed mechanism can reproduce accurate fluid and granular profiles. There are slight discrepancies in initial time instances, which may be attributed to instantaneous gate movement. However, upward/ downward gate removal (Fu and Jin, 2016) produces different results. A comprehensive evaluation of these issues is beyond the scope of the current study. It should also be noted that effective porosity is not applicable for numerical models having lengthscales smaller than the proper REV. In the present study, the minimum average interparticle spacing is assumed to be approximately
Fig. 12. Drag force distribution at 0.20 s for fluid particles (circle) and granular particles (triangle).
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
Fig. 13. Temporal variation of bed erosion.
6d, d being the granular particle diameter. Finer resolutions for continuum models may produce lesser numerical error, though they may not necessarily provide a physically accurate result (compared to their coarser counterparts). The initial effective porosity distribution for the existing/ proposed numerical model with coarse/ identical resolutions is shown in Fig. 11. Zero on the z-axis denotes the position of the granular media. The diffused effect is prominently visible for the coarse resolution for the erstwhile model. The effect of granular media quickly disappears outside the porous domain for either of the resolutions. The overestimation of effective porosity for fluid particles inside the granular media is reduced substantially for the proposed variant. For a superimposed fluid particle layer (at zeroth level), effective porosity does not vary considerably. However, zero distance is quite unlikely to happen as both modules are susceptible to rapid movement according to their corresponding governing equations. The distribution (in both present/ previous model) is heavily dependent on smoothing length and kernel choice. The larger
smoothing length (for coarse particle distribution) produces a diffused porosity distribution. A different approximation of kernel gradient (e.g. derivative approximation through MPS) may help for such cases. However, comparing various approximations of kernel gradients is beyond the scope of the current study. The drag force distribution at 0.25 s has been shown in Fig. 12 for different information mechanisms. Fluid velocity decreases suddenly at the granular surface for the current improved framework, unlike the previous case. Thus exposed granular particles are more prone to dislodgement due to the imposition of concentrated drag force instead of diffused distribution. The drag force attains its maximum value just outside the granular domain for the erstwhile framework. This defect is amplified by the combined errors in the calculation of effective porosity and particle number (np ). Calculation of relative velocity also influences the spread of drag forces far outside the granular media. Current modification produces a comparatively concentrated distribution of effective porosity and relative velocity. Thus, the drag force reduces rapidly outside the granular domain. The peak value
Fig. 14. Temporal variation of bed erosion for Scenario-III with identical (IIIa) and coarse (IIIb) sand distribution.
G. Pahar and A. Dhar
of drag force just coincides with the granular interface. The proposed model is essentially applicable for bed load transport. Some particles may get dislodged and move higher into fluid stream with large velocity. However, they do not undergo phase change or get diffused into fluid particles. The temporal variation of bed level is shown in Fig. 13. The bed level variation for third scenario are compared for existing and improved models with both identical and coarse resolutions. The patterns over 2 s have been depicted in Fig. 14. Even though the overall pattern for either of the cases remains similar, a few particles can be seen moving separately in the improved framework. It can be attributed to sharp velocity distribution of the fluid nearby granular interface. For the zeroth order consistent cases, the granular slope becomes stable within a short period (for Scenario-I, II). For the third scenario, similar trend can be observed for the improved/ previous models having both identical/ different resolutions. It may be noted that initially 5 horizontal layers of particle signify the granular domain (for Scenario-IIIb). Fig. 14 highlights the capability of the proposed model across different resolutions. The dislodgement of particles at the downstream of the gate can be related to dune migration. The proposed model is more suitable for such scenarios, whereas the previous model forms an apparent protective overlay of artificial drag force just above the granular interface. Figs. 9 and 10 show that the current variant can reproduce fluid and bed profiles satisfactorily. The current mechanism may slightly increase the computational time as an extra step involving gradient/divergence calculation has to be incorporated. Thus the corrective matrices need to be inverted at every location. However, searching of particles is required irrespective of the mechanisms. Overall the proposed modification consumes ≈ 4% of its predecessor. 4. Conclusion The current study demonstrates a significant improvement over existing coupled Lagrangian models (SPH-DEM, SPH-SPH) regarding applicability in lower/ identical resolutions. Inclusion of gradients reduces the diffused distribution of effective porosity. The information transfer mechanism has been tested for three cases of rapid sediment transport induced by instantaneous dam-break flows. The smeared effect of resistive forces arising out of granular media has been reduced substantially. Consequently, granular particles tend to travel longer distance while interacting with the rapidly moving fluid front. Effective porosity is calculated efficiently for the loosely disintegrated fluid-granular interface. The effect of artificial drag-force has been severely curtailed, though it cannot be eliminated completely. References Anderson, T.B., Jackson, R., 1967. Fluid mechanical description of fluidized beds. equations of motion. Ind. Eng. Chem. Fund. 6 (4), 527–539.
Advances in Water Resources 131 (2019) 103366 Bonet, J., Lok, T.-S., 1999. Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations. Comput. Methods Appl. Mech. Eng. 180 (1–2), 97–115. Boyce, C.M., Holland, D.J., Scott, S.A., Dennis, J.S., 2014. Novel fluid grid and voidage calculation techniques for a discrete element model of a 3d cylindrical fluidized bed. Comput. Chem. Eng. 65, 18–27. Canelas, R.B., Crespo, A.J., Domínguez, J.M., Ferreira, R.M., Gómez-Gesteira, M., 2016. SPH-DCDEM model for arbitrary geometries in free surface solidfluid flows. Comput. Phys. Commun. 202, 131–140. Di Felice, R., 1994. The voidage function for fluid-particle interaction systems. Int. J. Multiph. Flow 20 (1), 153–159. https://doi.org/10.1016/0301-9322(94)90011-6. Fu, L., Jin, Y.-c., 2016. Improved multiphase Lagrangian method for simulating sediment transport in dam-break flows. J. Hydraul. Eng. 142 (6), 1–13. Kafui, K., Thornton, C., Adams, M., 2002. Discrete particle-continuum fluid modelling of gassolid fluidised beds. Chem. Eng. Sci. 57 (13), 2395–2410. Kazemi, E., Nichols, A., Tait, S., Shao, S., 2017. SPH modelling of depth-limited turbulent open channel flows over rough boundaries. Int. J. Numer. Methods Fluids 83 (1), 3–27. Komoroczi, A., Abe, S., Urai, J.L., 2013. Meshless numerical modeling of brittle–viscous deformation: first results on boudinage and hydrofracturing using a coupling of discrete element method (DEM) and smoothed particle hydrodynamics (SPH). Comput. Geosci. 17 (2), 373–390. Lee, E.-S., Moulinec, C., Xu, R., Violeau, D., Laurence, D., Stansby, P., 2008. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. J. Comput. Phys. 227 (18), 8417–8436. https://doi.org/10.1016/j.jcp.2008.06.005. Markauskas, D., Kruggel-Emden, H., Sivanesapillai, R., Steeb, H., 2017. Comparative study on mesh-based and mesh-less coupled CFD-DEM methods to model particle-laden flow. Powder Technol. 305 (Supplement C), 78–88. Pahar, G., Dhar, A., 2017. Coupled incompressible Smoothed Particle Hydrodynamics model for continuum-based modelling of sediment transport. Adv. Water Resour. 102, 84–98. Pahar, G., Dhar, A., 2018. On force consideration in coupled isph framework for sediment transport in presence of free-surface flow. Environ. Fluid Mech. 18 (3), 555–557. Peng, Z., Moghtaderi, B., Doroodchi, E., 2016. A modified direct method for void fraction calculation in CFD-DEM simulations. Ad. Powder Technol. 27 (1), 19–32. Ren, B., Jin, Z., Gao, R., xue Wang, Y., lin Xu, Z., 2014. SPH-DEM modeling of the hydraulic stability of 2d blocks on a slope. J.Waterway Port Coast. Ocean Eng. 140 (6), 04014022. Robinson, M., Ramaioli, M., Luding, S., 2014. Fluidparticle flow simulations using two-way-coupled mesoscale SPH-DEM and validation. Int. J. Multiph. Flow 59 (Supplement C), 121–134. Shan, T., Zhao, J., 2014. A coupled CFD-DEM analysis of granular flow impacting on a water reservoir. Acta Mech. 225 (8), 2449–2470. Spinewine, B., 2005. Two-Layer Flow Behaviour and the Effects of Granular Dilatancy in Dam-Break Induced Sheet-Flow. Université catholique de Louvain Ph.D. thesis. Sun, R., Xiao, H., 2016. Sedifoam: a general-purpose, open-source CFD-DEM solver for particle-laden flow with emphasis on sediment transport. Comput. Geosci. 89, 207–219. Sun, X., Sakai, M., Yamada, Y., 2013. Three-dimensional simulation of a solidliquid flow by the DEM-SPH method. J. Comput. Phys. 248 (Supplement C), 147–176. Tan, H., Chen, S., 2017. A hybrid DEM-SPH model for deformable landslide and its generated surge waves. Adv. Water Resour.. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389–396. https://doi.org/10.1007/BF02123482. Wu, C., Zhan, J., Li, Y., Lam, K., Berrouk, A., 2009. Accurate void fraction calculation for three-dimensional discrete particle model on unstructured mesh. Chem. Eng. Sci. 64 (6), 1260–1266. Zhao, J., Shan, T., 2013. Coupled CFD DEM simulation of fl uid particle interaction in geomechanics. Powder Technol. 239, 248–258. Zhao, Y., Xu, L., Zheng, J., 2017. CFD-DEM simulation of tube erosion in a fluidized bed. AIChE J. 63 (2), 418–437.
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Short communication
On information coupling in hybrid ISPH framework for fluidized granular systems Gourabananda Pahar a,∗, Anirban Dhar b a b
Assistant Professor, Department of Civil Engineering, Indian Institute of Technology Kanpur, Kanpur UP 208016, India Associate Professor, Department of Civil Engineering, Indian Institute of Technology Kharagpur, Kharagpur WB 721302, India
a r t i c l e
i n f o
a b s t r a c t
Keywords: ISPH Sediment transport Information coupling
An improved information transfer mechanism is developed to couple fluid and granular modules for Incompressible Smoothed Particle Hydrodynamics modelling of sediment transport. The interaction force pair necessitates effective projection of information among parallel continuum based modules. The particle position in either module is susceptible to considerable variation due to their Lagrangian nature. An attempt has been made to minimize the effective thickness of the diffused interface between the fluid and granular particles. The proposed model reduces artificial drag force resulting from diffused computation of effective porosity. The revised framework outperforms in coarse resolution compared to its existing counterpart.
1. Introduction
the fraction of the volume associated with each fluid cell. Traditionally CFD-DEM assumes the full volume of the solid particle stays within the corresponding grid. Thus, all solid particles present within a fluid cell have equal contribution through their representative volume (irrespective of their positions). Peng et al. (2016) have introduced a Particle Meshing Method to take care of solid particle volume present within an Eulerian grid. In a Lagrangian method like SPH, the fluid cell volume can be mapped to the idea of influence domain. The contribution of particles has been generally taken through an uncorrected kernel value employed on them (Sun et al., 2013). The contribution of the particle decreases monotonically from the concerned location leading to an efficient representation. Typically, CFD-DEM is an ensemble of coupled Eulerian (CFD) and discrete (DEM) methods. The effective porosity and other coupling parameters are directly (or crudely) interpolated at the center of the fluid cell (Sun and Xiao, 2016; Zhao et al., 2017). Thus the resolution of fluid cells inevitably becomes coarser compared to the discrete granular module. In coupled SPH-DEM (Robinson et al., 2014; Markauskas et al., 2017; Tan and Chen, 2017) systems, particles in both parallel modules are moving separately. Thus the information needs to be suitably interchanged between modules. The drag force computation for any particle requires relative velocity at a common location, which is computed generally through a Shepard corrected kernel (Robinson et al., 2014). The pressure gradient and viscous term acting on the solid particle can be evaluated directly for a coupled CFD-DEM method from the corresponding fluid cell. Coupled SPH-DEM (Markauskas et al., 2017; Tan and Chen, 2017) generally avoids the pressure coupling by directly incorporating buoyancy force into account. It also removes the
Lagrangian particle methods are generally well-equipped to capture the interaction between rapid flow and fluidized granular bed (e.g., landslide-generated waves, sediment transport). Coupled flow models utilize various algorithms for interchanging information between submodules (.i.e. fluid and granular media). However, sharp discontinuity is often flattened out by these interpolation algorithms. This work aims to reduce the adverse impact of information loss by improving interpolation algorithms between different components of the coupled flow framework. The continuum conceptualization of effective porosity does not remain valid beyond a certain minimal representative elemental volume (REV). The volume-averaged flow equations lose their validity beyond this threshold length-scale corresponding to minimum REV. The distribution of effective porosity across interface should not depend on the resolution adopted for a particular simulation. Several researchers have formulated various elaborate algorithms related to Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) to achieve an actual distribution of effective porosity. Wu et al. (2009) have developed an efficient analytical method based on particle volume shared by various unstructured grids for coupled discrete particle methods. Generally, CFD-DEM directly calculates voidage function through the number of particles present in a particular fluid cell (Zhao and Shan, 2013; Shan and Zhao, 2014). Recent studies have suggested creating a timedependent mesh/ analytical solution (Boyce et al., 2014) based on current co-ordinates of granular particles at every time step to calculate
∗
Corresponding author. E-mail addresses: [email protected] (G. Pahar), [email protected] (A. Dhar).
https://doi.org/10.1016/j.advwatres.2019.06.011 Received 4 July 2018; Received in revised form 1 May 2019; Accepted 24 June 2019 Available online 26 June 2019 0309-1708/© 2019 Elsevier Ltd. All rights reserved.
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
likely to satisfy as particles tend to travel separately in their respective domains. Unlike the existing SPH-DEM formulations, the present framework utilizes initial identical (and coarser) resolution for the granular module. The improved technique utilizes a first-order consistent method to take care of rapid yet effective interpolation of information between parallel modules. Effect of artificial drag forces resulting in the previous model has been thoroughly minimized by the proposed variant. The modified framework is compared with existing experimental data along with numerical result simulated by the erstwhile formulation. 2. Information transfer between modules Anderson and Jackson (1967) had developed the averaged equations for fluid flow over deformable granular bed. Considering continuum approximation for the granular phase, the governing equations of continuity and momentum for granular and fluid phases can be written as (Pahar and Dhar, 2018), Fig. 1. Coupled modular structure.
possibility of enforcing the erroneous pressure gradient resulting from spurious numerical fluctuations in weakly compressible SPH. Generally, locally averaged Navier-Stokes equation requires a comparatively coarser resolution (compared to the granular module) of fluid phase to obtain sufficiently smooth porosity field (Sun et al., 2013; Tan and Chen, 2017). There are also instances of coupled SPH-DEM frameworks (Komoroczi et al., 2013; Ren et al., 2014; Canelas et al., 2016), where DEM resolution is coarser compared to SPH. However, they are not suitable for simulation of sediment-transport mechanism. A hybrid continuum based ISPH-ISPH framework has been developed in a previous study (Pahar and Dhar, 2017) for incompressible modelling of sediment transport. The framework shares similarities with existing SPH-DEM models though the granular media is modelled using an incompressible continuum conceptualization (Fig. 1) through rheological approximation. The modules interact with each other through an interaction force pair (Kafui et al., 2002). The force pair consists of fluid pressure-gradient, the divergence of viscous stress tensor and drag force. Each of these components would require extensive information from the respective parallel module. The fluid pressure is inherently smooth due to divergence-free nature of truly incompressible SPH (Lee et al., 2008). Thus fluid pressure gradient is incorporated into the interaction force pair without the additional assumption of hydrostatic pressure distribution (buoyancy). Different shades of blue in the continuum based fluid phase in Fig. 1 represent flow through porous media (lighter) and outside it (darker). The interface between two phases should be clearly delineated to distinguish the effect of resistive force vector originating from solid skeletal structures. The modules exist in parallel domains although their coordinate systems coexist. Interaction forces are calculated based on the projection of one module to another one. Prior work utilized a zeroth order smoothing algorithm to transmit information between modules. Though It ensures a smoothed distribution of different properties, the rapid discontinuity is flattened out leading to loss of information. For example, drag force varies non-linearly with corresponding effective porosity. It vanishes for zones having unit effective porosity. A large diffused interface between porous and non-porous zone would exert an extra artificial drag force on the fluid module. Thus it may adversely affect the computation of erosive processes. In the current study, an improved information transmission mechanism between coupled modules has been proposed. The projection of information becomes tricky for these separated yet coupled modules as both sets of particles are highly susceptible to irregular distribution. Gradient/ divergence(s) of properties are directly considered for projection in the higher order interpolation. SPH requires information regarding volume associated with the center particle to calculate the gradient/ divergence at the concerned location. However, these criteria are un-
∇ ⋅ 𝐮𝑠 = 0
(1a)
∇ ⋅ (𝜂𝐮𝑓 ) = 0
(1b)
𝐅𝑓 𝑠 𝐷𝐮𝑠 1 1 =− ∇𝑃𝑠 + 𝐠 + ∇ ⋅ 𝛕𝐬 + 𝐷𝑡 𝜙𝜌𝑠 𝜙𝜌𝑠 𝜙𝜌𝑠
(1c)
𝐷(𝜂𝐮𝑓 ) 𝐷𝑡
=−
𝐅 𝜂 1 ∇𝑃𝑓 + 𝜂𝐠 + 𝜈∇2 (𝜂𝐮𝑓 ) − 𝐷 + ∇ ⋅ 𝛕̄ 𝜌𝑓 𝜌𝑓 𝜌𝑓
(1d)
where 𝐮𝑠 = granular velocity, 𝐮𝑓 = fluid velocity, 𝜙 = solid fraction, 𝑃𝑠 = granular pressure, 𝐠 = acceleration due to gravity, 𝜌𝑠 = granular density, 𝐅𝑓 𝑠 = interaction force pair per unit volume, 𝜂 = effective porosity, 𝑃𝑓 = ̄ fluid pressure, 𝐅𝐷 = drag force per unit volume, 𝜈 = fluid viscosity, 𝜌𝛕 = 𝑓
sub-particle-scale tensor for modelling small-scale flow information. The interaction force can be represented as Kafui et al. (2002), 𝐅𝑓 𝑠 = −(1 − 𝜂)∇𝑃𝑓 + (1 − 𝜂)∇ ⋅ 𝛕 + 𝐅𝐷
(2)
Drag force (FD ) is determined through the following expression following Di Felice (1994). 𝐅𝐷 =
1 𝑛 𝐶 𝜌 𝜋𝑑 2 𝜂 (2−𝜉) (𝐮𝑓 − 𝐮𝑠 )|𝐮𝑓 − 𝐮𝑠 | 8 𝑝 𝑑𝑖 𝑓
(3)
where 𝑛𝑝 = number of granular particles present in a unit volume, 𝜉 = corrective porosity function. The drag force component is a non-linear function of effective porosity. In a continuum approach, np can be calculated as the volume of solid (1 − 𝜂) divided by the volume of a single particle ( 𝜋6 𝑑 3 ). Generally, bed-roughness is treated through drag-based formulation (Kazemi et al., 2017) near the bed-zone. The bed is considered as a rigid boundary for these models. The proposed model does not consider granular media as a solid boundary. Rather fluid outside and inside porous media is solved through a volume-averaged filtered Navier-Stokes equation. The momentum equation contains drag force which incorporates the information regarding particle diameter and relative velocity between fluid and granular particles. A diffused interface will ensure a smudged distribution of np , which in turn introduces artificial drag force even when the granular particles are not immediate neighbours of fluid particles. The information transfer mechanism between modules is described in detail in the following sections. 2.1. Information transfer: granular to fluid module The terms containing fluid pressure gradient and divergence of viscous stress are implicitly embedded into the fluid module. Thus, the informations required from the granular module for calculating drag force are the effective porosity and granular velocity at each fluid particle. According to Taylor Series, any quantity 𝜓 can be written as, | 𝜓𝑏 = 𝜓𝑎 + 𝐫𝑎𝑏 ⋅ ∇𝜓 | + ⋯ (4) |𝑎
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
where 𝐫𝑎𝑏 = relative position vector between particles a and b, 𝐫𝑎 = position vector of particle a. In discrete form the aforementioned equation can be written as, 𝜓𝑏 = 𝜓𝑎 + 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑ (𝜓𝑐 − 𝜓𝑎 )𝑣𝑐 ∇𝑤𝑎𝑐
(5)
𝑐=1
𝑐=1
Eq. (6) is valid for only a pair of particles inside the influence region. Similar sets of equations can be produced by sweeping over all the particles present inside support domain. Applying a kernel throughout the particles inside the support domain, the equation can be modified as, 𝑁𝑡 ∑
[(
𝑏=1 𝑁𝑡
=
∑ 𝑏=1
1 − 𝐫𝑎𝑏 ⋅ [(
| 𝜙𝑏 = 𝜙𝑎 + 𝐫𝑎𝑏 ⋅ ∇𝜙| |𝑎 = 𝜙𝑎 + 𝐫𝑎𝑏 ⋅
𝑐=1
𝑣𝑐 = volume associated with particle c present in the support domain of particle a, 𝑤𝑎𝑐 = kernel value at any particle c with respect to center particle a. The total number of particles present in a compact support domain (of a single module) is denoted by Nt . The center particle position is represented by a. However, the information regarding center particle position and volume associated with it remains unknown. All the particles present in the support domain except the concerned center particle remain in an identical module. The kernel gradient possesses zero value at the center location. Thus, Eq. (5) can be written as, ( ) ( ) 𝑁𝑡 𝑁𝑡 ∑ ∑ 𝜓𝑎 1 − 𝐫𝑎𝑏 ⋅ 𝑣𝑐 ∇𝑤𝑎𝑐 = 𝜓𝑏 − 𝐫𝑎𝑏 ⋅ 𝜓𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 (6)
𝜓𝑎
the domain, the solid fraction at central location can be written as,
𝑁𝑡 ∑ 𝑐=1
𝜓𝑏 − 𝐫𝑎𝑏 ⋅
)
]
𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏 𝑁𝑡
∑ 𝑐=1
)
(7)
2.1.1. Effective porosity calculation In CFD-DEM, effective porosity of a fluid cell is generally determined based on the total volume of solid particles present in the respective cell. To obtain the same in the coupled SPH-SPH model, a concept of ghost particles is introduced (Fig. 2). Fig. 2a depicts the representative particle positions for both fluid (blue) and granular (yellow) modules. The presence of granular particles in the support domain has been shown in Fig. 2b. A discontinuity of solid fraction exists exactly at the interface. The completeness of the corrupted support domain is maintained using a few virtual (grey) particles. The value of solid fraction associated with these particles is zero. The gradient (and higher order derivatives) of solid fraction is also considered to be zero. The location and associated volume of these ghost particles remain unknown. The representative volume associated with virtual granular particles need not be identical to that of actual granular particle. Outside the porous media, For the real particles present inside
Fig. 2. Effective porosity calculation.
𝑐=1 𝑁𝑡
∑ 𝑐=1
(𝜙𝑐 − 𝜙𝑎 )𝑣𝑐 ∇𝑤𝑎𝑐
(8)
𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐
Nt represents the total number of particles present in the compact support domain. Thus Nt becomes equal to 𝑁𝑎𝑐𝑡 + 𝑁𝑣𝑖𝑟 , where Nact and Nvir are respective numbers of actual and virtual particles inside the support domain. −𝜙𝑎 can be discarded as the influence domain has compact support (zeroth order consistency). For the virtual particles assumed inside the domain, the solid fraction at central location can be written as, | 𝜙𝑎 = 𝜙𝑏 + 𝐫𝑎𝑏 ⋅ ∇𝜙| |𝑏
(9)
Eq. (9) utilizes the information regarding the solid fraction and gradient of solid fraction (which are zero as discontinuity prevails outside granular domain) at the virtual particles. Accordingly, these particles have zero contribution at central location. The ghost particles also en∑𝑁 sures 𝑏 𝑡 𝑤𝑎𝑏 𝑣𝑏 to become unity for compact support domain. Thus solid fraction at any location can be written as, 𝜙𝑎 =
]
𝜓𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏
= 𝜙𝑎 + 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑
𝑁𝑎𝑐𝑡
∑
[(
𝑏=1
𝜙𝑏 − 𝐫𝑎𝑏 ⋅
𝑁𝑡 ∑ 𝑐=1
)
]
𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏
) ] ∑ [( | 𝜙𝑏 + 𝐫𝑎𝑏 ⋅ ∇𝜙| 𝑣𝑏 𝑤𝑎𝑏 |𝑏 𝑏=1 [( ) ] 𝑁𝑎𝑐𝑡 𝑁𝑎𝑐𝑡 ∑ ∑ 𝜙𝑏 − 𝐫𝑎𝑏 ⋅ 𝜙𝑐 𝑣𝑐 ∇𝑤𝑎𝑐 𝑣𝑏 𝑤𝑎𝑏 = 𝑁𝑣𝑖𝑟
+
𝑏=1
(10)
𝑐=1
For the virtual granular particles, solid fraction and gradient of solid fraction are zero. It may appear that zero solid-fraction and its gradient for virtual granular particles is rather conflicting. However, the calculation of gradient depends upon choice of smoothing length and representative volume associated with virtual particles. Considering very small smoothing length along with numerous virtual granular particles represent physically consistent result. Their number, location and associated
Fig. 3. Effective porosity distribution.
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
Fig. 4. Particle distribution.
Fig. 7. Information transfer: from fluid to granular media.
Fig. 5. Effective porosity distribution: (a) Direct averaging, (b) Direct averaging with edge trimming, (c) Zeroth order, (d) First Order, (e) actual.
Fig. 6. Error in calculation of effective porosity distribution: (a) Direct averaging, (b) Direct averaging with edge trimming, (c) Zeroth order, (d) First Order.
volumes do not feature in the equation. Thus, effective porosity at any point can be written as, 𝜂𝑎 = 1 − 𝜙𝑎
(11)
2.1.2. Velocity interpolation Granular velocity also needs to be interpolated on the fluid particle position for relative velocity calculation. Including the virtual particles in the simulation would require knowledge of their respective velocity, which is likely to be non-zero (unlike the previous case). Thus, the virtual particles cannot be directly included. The kernel gradient has to be corrected for corrupted support domain. The interpolated property 𝜓
Fig. 8. Initial conditions.
can be written as (Eq. (7)), ) ] ∑𝑁𝑎𝑐𝑡 [( ∑𝑁𝑎𝑐𝑡 𝜓𝑐 𝑣𝑐 ∇𝑤̃ 𝑎𝑐 𝑣𝑏 𝑤̃ 𝑎𝑏 𝜓𝑏 − 𝐫𝑎𝑏 ⋅ 𝑐=1 𝑏=1 𝜓𝑎 = ∑ [( ) ] ∑𝑁𝑎𝑐𝑡 𝑁𝑎𝑐𝑡 𝑣𝑐 ∇𝑤̃ 𝑎𝑐 𝑣𝑏 𝑤̃ 𝑎𝑏 1 − 𝐫𝑎𝑏 ⋅ 𝑐=1 𝑏=1
(12)
𝑤̃ and ∇𝑤̃ represent Shepard corrected kernel and corrected (not Shepard corrected) kernel gradient (Bonet and Lok, 1999). The kernel gradients do not need to be zeroth order consistent, as the difference
G. Pahar and A. Dhar
(𝜓𝑐 − 𝜓𝑎 ) is incorporated for computation of gradients. Zeroth and firstorder consistent kernel gradients will contain volume and location of the unknown central particle. Utilizing first order consistent kernels will implicitly interpolate the information of central particle. Representative distributions of effective porosity for different methods have been shown in Fig. 3. Fig. 3a depicts actual distribution of solid particles and corresponding effective porosity associated with fluid particles. The porosity can be calculated directly by averaging the void space in a unit volume (Fig. 3b). However, it leads to a diffused interface having width 𝜍, where 𝜍 denotes the diameter of the circular support domain. The partition of volume associated with fluid particles inside the support domain affects the calculation of effective porosity. The inverse distance weighting property of the kernel reduces the diffused nature of the interface (Zeroth order). Including the kernel gradients into the Eq. (7) would lead to a more rapid transition of effective porosity. Ideally, an abrupt change of effective porosity should exist between these zones. However, in a continuum approach, it is difficult to calculate the actual distribution of discontinuous effective porosity. Effective porosity cannot be calculated by utilizing the nearest neighbour approach (with minimal average in-
Advances in Water Resources 131 (2019) 103366
terparticle distance) due to its inherent continuity and differentiability constraints. To validate the proposed algorithm a representative coupled model has been shown. Distributions of fluid and granular particles have been depicted in Fig. 4. A regular particle distribution has been adopted for both modules with identical inter-particle distance. The topmost granular particle layer superimposes with the fluid particle layer. The distributions of effective porosity obtained from direct averaging (with and without volume trimming), kernel and kernel gradients have been shown in Fig. 5a–d. The exact distribution of effective porosity for the fluid particles has been depicted in Fig. 5e. The effective porosity of fluid particle layer coinciding with topmost granular layer is considered to be (1 − 𝜙2 ). Direct averaging considers identical contributions of particles present in a support domain. Each particle is considered to be associated with an identical circular volume (Fig. 5a). The effective porosity is calculated from the ratio of the total volume represented by the solid particles and the volume of the support domain (following traditional CFD-DEM methods). It underestimates the effective porosity for the interior portion of the granular domain. Total volume carried by every
Fig. 9. Effective porosity distribution at, (a) t = 0.25 s, (b) t = 0.50 s, (c) t = 0.75 s, (d) t = 1.00 s, (e) t = 1.25 s, (d) t = 1.50 s .
G. Pahar and A. Dhar
particle is included irrespective of its position. Total volume associated with the particles present in the vicinity of the boundary of the support domain may not remain within the support domain. Edges are trimmed for the second case (Fig. 5b). However, the width of the interface remains just as large as the previous case. The circular volume representation becomes essentially invalid for the current particle distribution. Each particle represents a square volume with the particle placed at the center of the square. A Voronoi diagram should be drawn at each timeinstance around the particles to get accurate volume representation as particle position is susceptible to large temporal variation. A straightforward kernel approximation is adopted in the third case (Fig. 5c). The corruption in the support domain is calculated through completeness of kernel. However, particle layers present above granular profile get heavily influenced. The current method provides suitable distribution of effective porosity leading to sharp interface (Fig. 5d). The error (Δ𝜂) in effective porosity calculation for each of these methods have been shown in Fig. 6, where Δ𝜂 = 𝜂exact − 𝜂calculated . The embedded kernel gradient influences the effective porosity through corresponding particle distance concerning the center particle. The kernel gradient peaks in the vicinity of the center particle (for both compact/ corrupted support domain). In SPH, a continuous kernel is transformed into discrete quantities for weight calculation. Thus a larger support domain would be a better choice regarding unity representation. However, the width of the diffused interface also increases with larger support domain. In the current study, 1.60Δx has been taken as smoothing length for a C2 kernel (Wendland, 1995). Maximum 6 layers of particles will be influenced by the proposed method for identical resolution of either module. The effect of kernel gradient vanishes for a full influence domain (the interior portion of granular media). The error in the calculation of effective porosity for a compact support domain is in order of 0.47%. Adopting higher order kernels (C4, C6) may further decrease the error in computation. However, being too concen-
Advances in Water Resources 131 (2019) 103366
trated around the central particle, they would require a larger support domain effectively increasing the width of the diffused interface. It should be also noted that the mechanism (for relative velocity calculation) is developed on the basis of zeroth order inconsistency ∑ 𝜕 𝑤̃ 𝑎𝑐 ( 𝑁 𝑐=1 𝜕𝑥 ≠ 0). Both Shepard corrected kernel and kernel gradient must include the center point and volume associated with it. In the current coupled modular framework, fluid and granular particles are present in the parallel domain while projecting information to another counterpart. Thus the volume of the projected center particle remains unknown. The volume of the center particle does not feature in the calculation of the proposed method. 2.2. Fluid to granular module Unlike the last case, the full interaction pair needs to be incorporated into the granular module. Thus the pressure gradient and divergence of stress have to be accurately calculated at every granular particle position. The procedure has been shown in the Fig. 7. The transfer process has been highlighted for a few isolated granular particles present inside the high-velocity fluid regime. Generally, the granular domain remains fully submerged leading to a compact support domain (in parallel modules). For erosion resulting from the highvelocity fluid stream, granular particles may become separated from the bed. The fluid property is calculated on the granular particles with zeroth order averaging through Shepard corrected kernel in the previous framework. Ultimately, the gradient/divergence of any property can be computed based on the support domain of the granular media in the final step. This procedure will be able to produce proper result inside the porous media; i.e., where each granular particle has a sufficient number of granular particles present in its support domain. However, for isolated granular particles, this may result in erroneous interpolation of information. The case shown in Fig. 7 contains 3–4 particles present
Fig. 10. Effective porosity distribution at, (a) t = 0.25 s, (b) t = 0.50 s, (c) t = 0.75 s, (d) t = 1.00 s for Scenarios-IIIa (identical resolution) and IIIb (coarse granular resolution).
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
divergence on each fluid particles. However, the procedure is completed within the first two steps of the aforementioned mechanism. 3. Results and discussions
Fig. 11. Initial distribution of effective porosity (𝜂) for Scenario-III.
in the influence zone each granular particle. Sparse irregular distribution of these particles (around the edge of the support domain) will result in a badly-scaled corrective matrix for kernel gradient modification (Bonet and Lok, 1999). Thus the kernel gradient may have to be taken as zero for these secluded particles. However, inside a highvelocity fluid stream, the pressure gradient is bound to be sufficiently high to affect the particle movement. Thus, a different procedure has been proposed in the current study. In the present work, the pressure is not directly interpolated into the granular particles. Instead, the fluid pressure gradient is calculated on each fluid particle. The pressure gradient is directly interpolated on the granular particles. The procedure requires additional computational overhead for searching and calculating pressure gradient/stress
The proposed mechanism is applied to three cases of dam-break induced erosion enclosed in a box. The initial condition is depicted in the Fig. 8. All three cases have sand as the bed material (Spinewine, 2005). The lock-gate is removed swiftly to study the temporal variation of bed. These cases are also simulated using the previous method and compared with the proposed variant. The resolution is deliberately kept coarse (initial interparticle spacing l0 = 0.0125 m) to assess the effectiveness of the model. All the simulation parameters (e.g. support domain radius, free-surface detection criteria, etc) have been kept identical for both variants. The third case is also replicated with separate resolutions for different modules (𝑙0 |𝑠𝑎𝑛𝑑 = 0.02𝑚, 𝑙0 |𝑤𝑎𝑡𝑒𝑟 = 0.0125𝑚) to study the effect of the developed algorithm. The smoothing length is taken as 1.60l0 . Calculation of effective porosity utilizes the smoothing length of the granular domain, whereas the transfer of fluid information to granular module considers smoothing length of the fluid module. In case of identical resolutions, the smoothing length for either of the cases remains unchanged. Simulated fluid profiles with the proposed framework along with experimental data have been shown in Fig. 9. The effective porosity serves as an indicator of the granular profile at different time instances. The experimental data of Spinewine (2005) has been shown by using circles (for water) and triangles (for granular profile). The progression of fluid fronts for the third scenario along with identical/ coarse resolutions (for granular media) has been depicted in Fig. 10 The coupled ISPH framework with the proposed mechanism can reproduce accurate fluid and granular profiles. There are slight discrepancies in initial time instances, which may be attributed to instantaneous gate movement. However, upward/ downward gate removal (Fu and Jin, 2016) produces different results. A comprehensive evaluation of these issues is beyond the scope of the current study. It should also be noted that effective porosity is not applicable for numerical models having lengthscales smaller than the proper REV. In the present study, the minimum average interparticle spacing is assumed to be approximately
Fig. 12. Drag force distribution at 0.20 s for fluid particles (circle) and granular particles (triangle).
G. Pahar and A. Dhar
Advances in Water Resources 131 (2019) 103366
Fig. 13. Temporal variation of bed erosion.
6d, d being the granular particle diameter. Finer resolutions for continuum models may produce lesser numerical error, though they may not necessarily provide a physically accurate result (compared to their coarser counterparts). The initial effective porosity distribution for the existing/ proposed numerical model with coarse/ identical resolutions is shown in Fig. 11. Zero on the z-axis denotes the position of the granular media. The diffused effect is prominently visible for the coarse resolution for the erstwhile model. The effect of granular media quickly disappears outside the porous domain for either of the resolutions. The overestimation of effective porosity for fluid particles inside the granular media is reduced substantially for the proposed variant. For a superimposed fluid particle layer (at zeroth level), effective porosity does not vary considerably. However, zero distance is quite unlikely to happen as both modules are susceptible to rapid movement according to their corresponding governing equations. The distribution (in both present/ previous model) is heavily dependent on smoothing length and kernel choice. The larger
smoothing length (for coarse particle distribution) produces a diffused porosity distribution. A different approximation of kernel gradient (e.g. derivative approximation through MPS) may help for such cases. However, comparing various approximations of kernel gradients is beyond the scope of the current study. The drag force distribution at 0.25 s has been shown in Fig. 12 for different information mechanisms. Fluid velocity decreases suddenly at the granular surface for the current improved framework, unlike the previous case. Thus exposed granular particles are more prone to dislodgement due to the imposition of concentrated drag force instead of diffused distribution. The drag force attains its maximum value just outside the granular domain for the erstwhile framework. This defect is amplified by the combined errors in the calculation of effective porosity and particle number (np ). Calculation of relative velocity also influences the spread of drag forces far outside the granular media. Current modification produces a comparatively concentrated distribution of effective porosity and relative velocity. Thus, the drag force reduces rapidly outside the granular domain. The peak value
Fig. 14. Temporal variation of bed erosion for Scenario-III with identical (IIIa) and coarse (IIIb) sand distribution.
G. Pahar and A. Dhar
of drag force just coincides with the granular interface. The proposed model is essentially applicable for bed load transport. Some particles may get dislodged and move higher into fluid stream with large velocity. However, they do not undergo phase change or get diffused into fluid particles. The temporal variation of bed level is shown in Fig. 13. The bed level variation for third scenario are compared for existing and improved models with both identical and coarse resolutions. The patterns over 2 s have been depicted in Fig. 14. Even though the overall pattern for either of the cases remains similar, a few particles can be seen moving separately in the improved framework. It can be attributed to sharp velocity distribution of the fluid nearby granular interface. For the zeroth order consistent cases, the granular slope becomes stable within a short period (for Scenario-I, II). For the third scenario, similar trend can be observed for the improved/ previous models having both identical/ different resolutions. It may be noted that initially 5 horizontal layers of particle signify the granular domain (for Scenario-IIIb). Fig. 14 highlights the capability of the proposed model across different resolutions. The dislodgement of particles at the downstream of the gate can be related to dune migration. The proposed model is more suitable for such scenarios, whereas the previous model forms an apparent protective overlay of artificial drag force just above the granular interface. Figs. 9 and 10 show that the current variant can reproduce fluid and bed profiles satisfactorily. The current mechanism may slightly increase the computational time as an extra step involving gradient/divergence calculation has to be incorporated. Thus the corrective matrices need to be inverted at every location. However, searching of particles is required irrespective of the mechanisms. Overall the proposed modification consumes ≈ 4% of its predecessor. 4. Conclusion The current study demonstrates a significant improvement over existing coupled Lagrangian models (SPH-DEM, SPH-SPH) regarding applicability in lower/ identical resolutions. Inclusion of gradients reduces the diffused distribution of effective porosity. The information transfer mechanism has been tested for three cases of rapid sediment transport induced by instantaneous dam-break flows. The smeared effect of resistive forces arising out of granular media has been reduced substantially. Consequently, granular particles tend to travel longer distance while interacting with the rapidly moving fluid front. Effective porosity is calculated efficiently for the loosely disintegrated fluid-granular interface. The effect of artificial drag-force has been severely curtailed, though it cannot be eliminated completely. References Anderson, T.B., Jackson, R., 1967. Fluid mechanical description of fluidized beds. equations of motion. Ind. Eng. Chem. Fund. 6 (4), 527–539.
Advances in Water Resources 131 (2019) 103366 Bonet, J., Lok, T.-S., 1999. Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations. Comput. Methods Appl. Mech. Eng. 180 (1–2), 97–115. Boyce, C.M., Holland, D.J., Scott, S.A., Dennis, J.S., 2014. Novel fluid grid and voidage calculation techniques for a discrete element model of a 3d cylindrical fluidized bed. Comput. Chem. Eng. 65, 18–27. Canelas, R.B., Crespo, A.J., Domínguez, J.M., Ferreira, R.M., Gómez-Gesteira, M., 2016. SPH-DCDEM model for arbitrary geometries in free surface solidfluid flows. Comput. Phys. Commun. 202, 131–140. Di Felice, R., 1994. The voidage function for fluid-particle interaction systems. Int. J. Multiph. Flow 20 (1), 153–159. https://doi.org/10.1016/0301-9322(94)90011-6. Fu, L., Jin, Y.-c., 2016. Improved multiphase Lagrangian method for simulating sediment transport in dam-break flows. J. Hydraul. Eng. 142 (6), 1–13. Kafui, K., Thornton, C., Adams, M., 2002. Discrete particle-continuum fluid modelling of gassolid fluidised beds. Chem. Eng. Sci. 57 (13), 2395–2410. Kazemi, E., Nichols, A., Tait, S., Shao, S., 2017. SPH modelling of depth-limited turbulent open channel flows over rough boundaries. Int. J. Numer. Methods Fluids 83 (1), 3–27. Komoroczi, A., Abe, S., Urai, J.L., 2013. Meshless numerical modeling of brittle–viscous deformation: first results on boudinage and hydrofracturing using a coupling of discrete element method (DEM) and smoothed particle hydrodynamics (SPH). Comput. Geosci. 17 (2), 373–390. Lee, E.-S., Moulinec, C., Xu, R., Violeau, D., Laurence, D., Stansby, P., 2008. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. J. Comput. Phys. 227 (18), 8417–8436. https://doi.org/10.1016/j.jcp.2008.06.005. Markauskas, D., Kruggel-Emden, H., Sivanesapillai, R., Steeb, H., 2017. Comparative study on mesh-based and mesh-less coupled CFD-DEM methods to model particle-laden flow. Powder Technol. 305 (Supplement C), 78–88. Pahar, G., Dhar, A., 2017. Coupled incompressible Smoothed Particle Hydrodynamics model for continuum-based modelling of sediment transport. Adv. Water Resour. 102, 84–98. Pahar, G., Dhar, A., 2018. On force consideration in coupled isph framework for sediment transport in presence of free-surface flow. Environ. Fluid Mech. 18 (3), 555–557. Peng, Z., Moghtaderi, B., Doroodchi, E., 2016. A modified direct method for void fraction calculation in CFD-DEM simulations. Ad. Powder Technol. 27 (1), 19–32. Ren, B., Jin, Z., Gao, R., xue Wang, Y., lin Xu, Z., 2014. SPH-DEM modeling of the hydraulic stability of 2d blocks on a slope. J.Waterway Port Coast. Ocean Eng. 140 (6), 04014022. Robinson, M., Ramaioli, M., Luding, S., 2014. Fluidparticle flow simulations using two-way-coupled mesoscale SPH-DEM and validation. Int. J. Multiph. Flow 59 (Supplement C), 121–134. Shan, T., Zhao, J., 2014. A coupled CFD-DEM analysis of granular flow impacting on a water reservoir. Acta Mech. 225 (8), 2449–2470. Spinewine, B., 2005. Two-Layer Flow Behaviour and the Effects of Granular Dilatancy in Dam-Break Induced Sheet-Flow. Université catholique de Louvain Ph.D. thesis. Sun, R., Xiao, H., 2016. Sedifoam: a general-purpose, open-source CFD-DEM solver for particle-laden flow with emphasis on sediment transport. Comput. Geosci. 89, 207–219. Sun, X., Sakai, M., Yamada, Y., 2013. Three-dimensional simulation of a solidliquid flow by the DEM-SPH method. J. Comput. Phys. 248 (Supplement C), 147–176. Tan, H., Chen, S., 2017. A hybrid DEM-SPH model for deformable landslide and its generated surge waves. Adv. Water Resour.. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv. Comput. Math. 4 (1), 389–396. https://doi.org/10.1007/BF02123482. Wu, C., Zhan, J., Li, Y., Lam, K., Berrouk, A., 2009. Accurate void fraction calculation for three-dimensional discrete particle model on unstructured mesh. Chem. Eng. Sci. 64 (6), 1260–1266. Zhao, J., Shan, T., 2013. Coupled CFD DEM simulation of fl uid particle interaction in geomechanics. Powder Technol. 239, 248–258. Zhao, Y., Xu, L., Zheng, J., 2017. CFD-DEM simulation of tube erosion in a fluidized bed. AIChE J. 63 (2), 418–437.