SPH-DEM Examples
Subaqueous granular collapse with a sphere
The coupling of a three-phase case can be set up using the SPH-DEM coupling. The submerged granular collapse shows the interaction of fluid and soil. Besides, a sphere is included to show its interaction with the other two phases. The SPH points used to represent fluid and soil overlap, unlike the DEM particles that are separated by their borders.
// MechSys
#include <mechsys/sph/Domain.h>
#include <iostream>
#include <fstream>
using std::cout;
using std::endl;
double DampF,DampS;
double tset = 0.1;
void UserFunctions(SPH::Domain &doma){
if (doma.Time<=tset){
#pragma omp parallel for schedule (static) num_threads(doma.Nproc)
for (size_t i=0; i<doma.Particles.Size(); i++){
if ( (doma.Particles[i]->IsFree) && (doma.Particles[i]->Material == 1 ) ){ doma.Particles[i]->a -= DampF * doma.Particles[i]->v;}
else if ( (doma.Particles[i]->IsFree) && (doma.Particles[i]->Material == 3 ) ){ doma.Particles[i]->a -= DampS * doma.Particles[i]->v;}
}
}
}
int main(int argc, char **argv) try
{
size_t Nproc = omp_get_max_threads();
if (argc>1) Nproc = atoi(argv[1]);
SPH::Domain dom;
dom.Nproc = Nproc;
dom.Dimension = 2;
dom.Scheme = 0;
double g = 9.81;
dom.Gravity = 0.0,-g,0.0;
dom.KernelType = 1;
dom.GradientType= 0;
dom.VisEq = 0;
dom.GeneralAfter = &UserFunctions;
dom.XSPH = 0.5;
dom.DeltaSPH = 0.15;
//Simulation time
double Tfinal = 2.0;
int numfiles = 50;
double dtoutput = Tfinal/numfiles;
char file_name[10] = "sub";
// SIZE AND CONDITIONS FOR THE DOMAIN
double L = 0.25; // Length of the tank
double H = 0.15; // Height of the tank
double Lc = 0.10; // Length of the granular column
double Hc = 0.10; // Height of the granular column
double n = 50; // number of points to discretise the domain
double dx = H/n;
double h = dx*1.3;
dom.InitialDist = dx;
double CFL = 0.05;
// PROPERTIES FOR FLUID
double rhoF = 1000.0; // Density of fluid
double Mu = 1.0e-3; // Dynamic viscosity
double DrhoF = 0.1*rhoF; // Delta density, allowed density variation for fluid.
double CsF = sqrt((g*L*rhoF)/DrhoF); // Speed of sound for fluid
double P0 = 0.0; // Background pressure
double Tf = (CFL*h/CsF); // Time-step for fluid
DampF = 0.2*CsF/h; // Damping coefficient for fluid
// PROPERTIES FOR SOIL
double rhoS = 1900.0; // Density of soil
double E = 10.0e6; // Young modulus
double Nu = 0.30; // Poisson's ratio
double Phi = 27.0; // Angle of shearing resistance
double Psi = 0.0; // Dilatancy angle
double c = 0.0; // Cohesion
double K = E/(3.0*(1.0-2.0*Nu)); // Bulk modulus
double G = E/(2.0*(1.0+Nu)); // Shear modulus
double ns = 0.4; // Porosity
double k = 1.0e-8; // Intrinsic permeability
double CsS = sqrt( (4.0*G)/(3.0*rhoS) + (K/rhoS) ); // Speed of sound for soil
double Ts = ((CFL*h)/CsS); // Time-step for the soil
DampS = 0.2*sqrt(E/(rhoS*h*h)); // Damping coefficient for soil
// PROPERTIES AND LOCATION FOR THE DEM SPHERE OR BOULDER
double Cx = 0.75*Lc; // Centre of sphere in y direction
double Cy = 0.75*Hc; // Centre of sphere in x direction
double Cz = 0.0; // Centre of sphere in x direction
double Rc = 0.05*H; // Sphere radius
double rhoB = 2200.0; // Density of sphere
// SOIL-FLUID INTERACTION
dom.SWIType = 2;
double timestep = std::min({Tf,Ts}); // Time-step Selected among fluid and soil time-steps
double steps = Tfinal/timestep;
//////////////////////////////////////// SPH PARTICLES CREATION ///////////////////////////////////////
dom.AddBoxLength(/*Tag*/1 ,/*Position*/OrthoSys::O, /*Length*/L , /*Height*/H , /*Width*/0.0, /*Radius*/dx/2.0 ,/*Density*/rhoF, /*Smoothing*/h, /*Packing*/1 , /*Rotation*/0 , /*Random*/false, /*Fixed*/false ); // FLUID
dom.AddBoxLength(/*Tag*/2 ,/*Position*/OrthoSys::O, /*Length*/Lc, /*Height*/Hc, /*Width*/0.0, /*Radius*/dx/2.0 ,/*Density*/rhoS, /*Smoothing*/h, /*Packing*/1 , /*Rotation*/0 , /*Random*/false, /*Fixed*/false ); // SOIL
for (size_t i=0; i<dom.Particles.Size(); i++){
double xb=dom.Particles[i]->x(0);
double yb=dom.Particles[i]->x(1);
if (dom.Particles[i]->ID==1){
// ASSIGNING MATERIAL PROPERTIES TO THE FLUID PARTICLES//
dom.Particles[i]->Material = 1;
dom.Particles[i]->Constitutive = 11;
dom.Particles[i]->Mu = Mu;
dom.Particles[i]->MuRef = Mu;
dom.Particles[i]->PresEq = 0;
dom.Particles[i]->P0 = P0;
dom.Particles[i]->Cs = CsF;
dom.Particles[i]->Alpha = 0.05;
dom.Particles[i]->Beta = 0.05;
dom.Particles[i]->TI = 0.3;
dom.Particles[i]->TIn = 4.0;
dom.Particles[i]->TIInitDist = dx;
dom.Particles[i]->Shepard = true;
}
else if (dom.Particles[i]->ID==2){
// ASSIGNING MATERIAL PROPERTIES TO THE SOIL PARTICLES//
dom.Particles[i]->Material = 3;
dom.Particles[i]->Constitutive = 31;
dom.Particles[i]->Fail = 3;
dom.Particles[i]->c = c;
dom.Particles[i]->phi = (Phi*M_PI)/180.0;
dom.Particles[i]->psi = (Psi*M_PI)/180.0;
dom.Particles[i]->K = K;
dom.Particles[i]->G = G;
dom.Particles[i]->SeepageType = 0;
dom.Particles[i]->k = k;
dom.Particles[i]->n = ns;
dom.Particles[i]->RhoF = 0.0;
dom.Particles[i]->Cs = CsS;
dom.Particles[i]->Alpha = 1.0;
dom.Particles[i]->Beta = 1.0;
dom.Particles[i]->TI = 0.5;
dom.Particles[i]->TIn = 2.55;
dom.Particles[i]->TIInitDist = dx;
dom.Particles[i]->Shepard = true;
}
// To delete some SPH particles and to put a DEM sphere into that empty space:
if ( sqrt( pow(xb-Cx,2.0) + pow(yb-Cy,2.0) )<(Rc+0.5*dx) ) dom.Particles[i]->ID = 0;
}
// To delete particles with an ID equals zero.
dom.DelParticles(0);
dom.DomMax = Vec3_t(L+2.0*dx,H+0.2*H,0.0);
dom.DomMin = Vec3_t(-L-2.0*dx,-2.0*dx,0.0);
//////////////////////////////////////// DEM PARTICLES CREATION ///////////////////////////////////////
dom.AddSegment(/*Tag*/ -1,/*PositionB*/Vec3_t(-dx,-dx,0.0),/*PositionT*/Vec3_t(-dx,H+0.1*H,0.0),/*Thickness*/dx,/*Density*/rhoS); // Left wall
dom.AddSegment(/*Tag*/ -2,/*PositionB*/Vec3_t(L+dx,-dx,0.0),/*PositionT*/Vec3_t(L+dx,H+0.1*H,0.0),/*Thickness*/dx,/*Density*/rhoS); // Right wall
dom.AddSegment(/*Tag*/ -3,/*PositionL*/Vec3_t(-dx,-dx,0.0),/*PositionR*/Vec3_t(L+dx,-dx,0.0),/*Thickness*/dx,/*Density*/rhoS); // Bottom wall
dom.AddSphere(/*Tag*/ -4, /*Centre*/Vec3_t(Cx,Cy,Cz),/*Radius*/Rc,/*Density*/rhoB); // Sphere
for (size_t j = 0;j<dom.DEMParticles.Size();j++){
dom.DEMParticles[j]->Ff = dom.DEMParticles[j]->Props.m*dom.Gravity;
dom.DEMParticles[j]->Props.Mu = tan((Phi*M_PI)/180.0); // friction coefficient in DEM segments
if ( dom.DEMParticles[j]->Tag == -4 ){continue;} // if sphere then avoid the command lines below
dom.DEMParticles[j]->FixVeloc(); // to hold the DEM segments so that they will not be affected by the gravity
dom.DEMParticles[j]->Props.eps = 0.0; // overlaping distance
for (size_t i=0; i<dom.Particles.Size(); i++){
double dist = DEM::Distance(dom.Particles[i]->x,*dom.DEMParticles[j]->Faces[0]);
if (dist<dom.DEMParticles[j]->Props.eps || dom.DEMParticles[j]->Props.eps == 0.0) dom.DEMParticles[j]->Props.eps = dist;
}
}
dom.DEMParticles[0]->Props.Mu = 0.0;
dom.DEMParticles[3]->Props.eps = 0.5*dx;
dom.DEMstiff = 0.002;
std::cout<<"CsF = " <<CsF<<std::endl;
std::cout<<"CsS = " <<CsS<<std::endl;
std::cout<<"Steps = " <<steps<<std::endl;
std::cout<<"Resolution Sphere-Fluid = "<<(2.0*Rc/dx)<<std::endl;
dom.Solve(/*tf*/Tfinal,/*dt*/timestep,/*dtOut*/dtoutput,file_name,numfiles);
return 0;
}
MECHSYS_CATCH
Initialisation
MechSys/SPH-DEM simulations are prepared with a C++ source code. The first step is to include the required MechSys libraries. In this case, the file domain.h contains the DEM domain definitions:
#include <mechsys/sph/domain.h>
Simulation script structure
The simulation script structure is:
#include ...
int main(int argc, char **argv) try
{
...
}
MECHSYS_CATCH
After including the required MechSys' libraries (and C++ standard libraries, if necessary), the main function is declared.
For MechSys, it is essential to add the try keyword just after the int main() and before the first opening brace {. Then, after the last closing brace, the MECHSYS_CATCH is added.
Declaring the SPH Domain
The following line returns an upper bound on the number of threads that could be used:
size_t Nproc = omp_get_max_threads();
If the first argument is introduced following the execution command, the number of threads will be reassigned:
if (argc>1) Nproc = atoi(argv[1]);
The line,
SPH::Domain dom;
declares the container dom of the class SPH::Domain. This class defines the SPH universe and has some functions to include SPH particles and DEM particles, some of which will be used in this tutorial.
The SPH model has many options to solve the equations. Thus, once the domain is named (e.g., dom), the functions and their options are selected using such an extension. Hence, the dimensionality of the problem can be selected using an integer, 2 or 3:
dom.Dimension = 2;
There are two explicit schemes to solve the temporal derivatives:
- 0 → Modified Verlet (Ghadimi, 2012)
- 1 → Leapfrog (Braune & Lewiner, 2013)
And it can be called as follow:
dom.Scheme = 0;
The gravity is a vector defined in the domain, which directions are x, y and z. A double variable was defined for the magnitude of gravity in the previous line, double g=9.81;.
dom.Gravity = 0.0,-g,0.0;
Many kernels can be found in the literature. Mechsys has three options (Korzani et al., 2017):
- 0 → Cubic kernel
- 1 → Quintic
- 2 → Quintic Spline
dom.KernelType = 1;
Remark: The option called Quintic is known as Wendland C2 or Wendland Quintic kernel function. Wendland kernel avoids pairing (or clumping) instability that occurs where there are more neighbouring SPH particles in the influence domain that a kernel function can accommodate (Dehnen & Aly, 2012; Bui & Nguyen, 2020).
Two options are available to discretise the gradient of the pressure or the stress tensor depending on the form of the density as a denominator (Monaghan, 2005):
- 0 → Squared density
- 1 → Multiplied density
dom.GradientType = 0;
If the system to be solved is a Newtonian fluid, there are four options to discretise the second-order derivative term:
- 0 → Morris et al. (1997).
- 1 → Shao & Lo (2003).
- 2 → Real viscosity for incompressible fluids.
- 3 → Real viscosity for compressible fluids (Takeda et al., 1994).
Then, it can be choose as:
dom.VisEq = 0;
There is a function in SPH that allows the user to program their own functions being included in the main code before or after computing the interaction forces. It can be done by using the following command, dom.GeneralAfter = &FunctionsName; or dom.GeneralBefore = &FunctionsName;. For instance:
dom.GeneralAfter = &UserFunctions;
The FunctionsName, in this case, is defined before the main structure of the code as a separated function. void type function from c++ is used for this purpose, such as shown below.
double DampF,DampS;
double tset = 0.1;
void UserFunctions(SPH::Domain &doma){
if (doma.Time<=tset){
#pragma omp parallel for schedule (static) num_threads(doma.Nproc)
for (size_t i=0; i<doma.Particles.Size(); i++){
if ( (doma.Particles[i]->IsFree) && (doma.Particles[i]->Material == 1 ) ){ doma.Particles[i]->a -= DampF * doma.Particles[i]->v;}
else if ( (doma.Particles[i]->IsFree) && (doma.Particles[i]->Material == 3 ) ){ doma.Particles[i]->a -= DampS * doma.Particles[i]->v;}
}
}
}
The user's function, in this case, will be used to set up another type of dissipative filter (artificial viscous damping) proposed to diminish the noise of soils stresses (Nguyen et al., 2017). However, a more forceful purpose is to cushion the sudden activation of gravity that ends up in particles falling simultaneously. Compared to many experiments, it might result in an artificial fact since the materials are already at rest. Thus, the activation of the filter damps the force with what particles shocks with others in the first moments of the simulation while the case reaches a steady-state or static condition. Also, it is possible to configure other functions such as
movable walls or gates, specifying their velocity or new positions at each time-step by using dom.GeneralBefore = &FunctionsName;. The configuration of the two parameters DampF or DampS can be found in Nguyen et al., (2017) and below. The activation is done by a brief period defined as double tset = 0.1;.
Also, it is possible to use some correction expressions to improve the solution given by SPH. XSPH is a widely used term to correct the position of the particles. This can be activated by calling the following function and assigning a value between 0 and 1. A typical value is 0.5:
dom.XSPH = 0.5;
The basic SPH structure produces a checkerboard pattern so that a filter to smooth the pressure field is necessary. Thus, Delta-SPH is a diffusive term added in the mass conservation equation that can be activated by calling the following function and assigning a value between 0 and 1, depending on the case.
dom.DeltaSPH = 0.15;
Then Tfinal defines the final time of the simulation and numfiles the number of files to print the solution of the system. dtoutput defines the time-step to print the files, but it is not related to the real time-step and stability condition. The line char file_name[10] = "sub"; sets up the mane of the output files, called sub in this case.
double Tfinal = 2.0;
int numfiles = 50;
double dtoutput = Tfinal/numfiles;
char file_name[10] = "sub";
Some geometrical parameters are defined to configure the dimensions of the problem. Thus, height H and length L of a box that will contain the SPH fluid particles are defined. Besides, height Hc and length Lc of a column that will contain the SPH soil particles are defined as well. The number of particles n, the spacing among them dx and the smoothing length double h = dx*1.3 are configured.
double H = 0.15; // Height of the tank
double L = 0.25; // Length of the tank
double Hc = 0.10; // Height of the granular column
double Lc = 0.10; // Length of the granular column
double n = 50; // number of points to discretise the domain
double dx = H/n;
double h = dx*1.3;
The value 1.3 is a scaling factor of the spread of the smoothing function, and it can be from 0.6 to 2.0; however, values between 1.1 and 1.3 are highly recommended. Also, it is necessary to assign the distance between points using the following function:
dom.InitialDist = dx;
The stability condition for both fluid and soil is defined by the Courant-Friedrichs-Lewy condition (CFL). Values lower than 0.1 will ensure the stability of the solution, depending on the study case.
double CFL = 0.05;
Then, some of the variables or parameter that can influence the time-step and the stability of the fluid are defined. A few parameters are need if the fluid is assumed as Newtonian: density rhoF, dynamic viscosity Mu, the allowed density variation of the fluid that is 10% in this case double DrhoF = 0.1*rhoF;, speed of the sound of the fluid CsF, a background pressure P0, the time-step Tf and the damping coefficient DampF for the fluid:
double rhoF = 1000.0; // Density of fluid
double Mu = 1.0e-3; // Dynamic viscosity
double DrhoF = 0.1*rhoF; // Delta density, allowed density variation for fluid.
double CsF = sqrt((g*L*rhoF)/DrhoF); // Speed of sound for fluid
double P0 = 0.0; // Background pressure
double Tf = (CFL*h/CsF); // Time-step for fluid
DampF = 0.2*CsF/h; // Damping coefficient for fluid
Remark: The time-step is defined as a function of the CFL. The commonly chosen sound speed for simulations differs from the physical value to make the solution finishable without affecting the results. Otherwise, the simulation might be endless either for fluid or soil. Morris (1997) shows some definitions regarding the speed of the sound and stability condition.
Then, some of the variables or parameter that can influence the time-step and the stability of the soil are defined. There a re two models for representing the soil, elpastoplatic and hypoplastic. The first option implemented in Mechsys is proposed by Bui et al. (2008) that has been implemented in many papers. It is an elastic perfectly-plastic model with a Drucker-Prager failure criterion, and accordingly, to the constitutive law, specific parameters are needed. Thus, the following parameters are setup: density of soil rhoS, cohesion c, angle of shearing resistance Phi, dilatancy angle Psi, Yong modulus E, Poisson's ratio Nu, bulk modulus K, shear modulus G, porosity ns, intrinsic permeability k, speed of the sound for soil CsS, the time-step Ts and the damping coefficient DampS for the soil:
double rhoS = 1900.0; // Density of soil
double c = 0.0; // Cohesion
double Phi = 27.0; // Angle of shearing resistance
double Psi = 0.0; // Dilatancy angle
double E = 10.0e6; // Young modulus
double Nu = 0.30; // Poisson's ratio
double K = E/(3.0*(1.0-2.0*Nu)); // Bulk modulus
double G = E/(2.0*(1.0+Nu)); // Shear modulus
double ns = 0.4; // Porosity
double k = 1.0e-8; // Intrinsic permeability
double CsS = sqrt( (4.0*G)/(3.0*rhoS) + (K/rhoS) ); // Speed of sound for soil
double Ts = ((CFL*h)/CsS); // Time-step for the soil
DampS = 0.2*sqrt(E/(rhoS*h*h)); // Damping coefficient for soil
If some soil particles whose diameter is much larger than the average, they can be represented by DEM spheres, called here "boulder". One DEM sphere is configured in this case. The centre of the DEM particle Cx, Cy, Cz, its radius Rc and density rhoB are needed.
double Cx = 0.75*Lc; // Centre of sphere in y direction
double Cy = 0.75*Hc; // Centre of sphere in x direction
double Cz = 0.0; // Centre of sphere in x direction
double Rc = 0.05*H; // Sphere radius
double rhoB = 2200.0; // Density of boulder
Since we have fluid and soil SPH particles interacting, the type of seepage force is needed to be defined. Mechsys has two options (0 and 1) defined for static conditions due to the use of the buoyancy concept. Option 2 has the seepage and the pore fluid pressure used for this case. And, option 3 inhibits the interaction between fluid and soil SPH particles. The seepage force is chosen by the command shown below.
- 0 → Seepage + buoyancy (Korzani et al., 2018)
- 1 → Seepage + surface erosion + buoyancy (Korzani et al., 2018)
- 2 → Seepage + pore fluid pressure (Bui et al., 2007)
- 3 → No fluid-soil interaction
dom.SWIType = 2;
There are two time-steps, one for fluid Tf and one for soil Ts, as shown above. The whole system's stability criteria are defined by the minimum time-step required for the SPH particles, either fluid or soil, through the following command double timestep = std::min({Tf,Ts});. The number of steps for the simulation can be computed as double steps = Tf/timestep;.
double timestep = std::min({Tf,Ts});
double steps = Tf/timestep;
Declaring the SPH Particles
A function that will insert the SPH particles into the domain is called. There are three functions that can be used for such a purpose that are AddBoxLength, AddSingleParticle and AddBoxNo, and the first one will be implemented in this example. The case programmed here is a multiphase flow, so that, two layers of SPH particles are configured with different size each. The function requires some parameters described inside the parentesis:
dom.AddBoxLength(/*Tag*/1 ,/*Position*/OrthoSys::O, /*Length*/L , /*Height*/H , /*Width*/0.0, /*Radius*/dx/2.0 ,/*Density*/rhoF, /*Smoothing*/h, /*Packing*/1 , /*Packing Rotation*/0 , /*Random*/false, /*Fixed*/false ); // FLUID
dom.AddBoxLength(/*Tag*/2 ,/*Position*/OrthoSys::O, /*Length*/Lc, /*Height*/Hc, /*Width*/0.0, /*Radius*/dx/2.0 ,/*Density*/rhoS, /*Smoothing*/h, /*Packing*/1 , /*Packing Rotation*/0 , /*Random*/false, /*Fixed*/false ); // SOIL
If the other two options are needed to insert SPH particles, their structures are given as follows:
dom.AddBoxNo( Tag, X, nx, ny, nz, r, Density, h, Packing, Packing rotation, random, Fixed);
Where nx, ny and nz are the numbers of SPH particles in each direction x, y and z, respectively.
To add a single particle, use:
dom.AddSingleParticle( Tag, X, Mass, Density, h, Fixed);
The first parameter is a "Tag" to distinguish different subgroups of particles. The vector X=Vec3_t(0.0, 0.0, 0.0) indicates the origin of the box with the SPH particles,which coordinates can be seen as Vec3_t(Left, Bottom, Front). Then, "Length", "Height", and "Width" must be defined by a double. The "Radius" is typically defined as dx/2. Then, two parameters are called, Density for each material and smoothing length, h. "Rotation": this option gets an integer value to define the direction in Hexagonal Close Packing. "Random": it is a boolean variable to add randomness in the position of each particle proportional to its radius. "Fixed": it is a boolean variable for defining the fixity of a particle.
Remark: rectangle and cube are the primary shapes discretised by the SPH points, in 2D and 3D, respectively. Then, the rectangle and cube can be reshaped by deleting some particles placed in unwanted zones using basic functions (e.g., the equation of the straight line, plane, circle, sphere, etc.). It can be done by changing the initial tag of the particles through the variable called dom.Particles[i]->ID = value. The particles with that new ID can be eliminated using a command from the code called dom.DelParticles(ID) as will be shown below.
After calling the function to create the SPH points, it is necessary to assign the parameters to each particle that define the type and behaviour of the material. This will be done by using a loop for, which structure is written as:
for (size_t i=0; i <dom.Particles.Size(); i++){
...
}
The coordinates can be re-assigned to facilitate the construction of the model by using conditionals. The re-assignation can be done via the following commands:
double xb=dom.Particles[i]->x(0);
double yb=dom.Particles[i]->x(1);
It is necessary to use the following structure to assign or define certain properties: dom.Particles[i]->parameter = value. However, if there are two types of materials a conditional based on the Tag can be used to assign the properties to each material as shown below:
if (dom.Particles[i]->ID==1){ \\ Fluid properties
...
}
else if (dom.Particles[i]->ID==2){ \\ Soil properties
...
}
Thus, the first parameter refers to the type of material, where there are three options that can be selected by using an integer:
- 1 → Fluid
- 2 → Solid
- 3 → Soil
dom.Particles[i]->Material = 1;
There are also some options to choose the constitutive model. There are five options for fluids, one for solids and two for soils:
- 11 → Linear (Newton, 1687).
- 12 → Bingham regularised (Papanastasiou, 1987).
- 13 → Bingham-Bilinear (Imran et al., 2001).
- 14 → Cross model (Shao & Lo 2003).
- 15 → Herschel-Bulkley (Ké & Turcotte, 1980).
- 2 → Elastic Perfect-Plastic with Von-Mises failure criteria.
- 31 → Elastic Perfect-Plastic with Drucker-Praguer failure criteria (Bui et al., 2008).
- 32 → Hypoplastic model (Peng et al., 2015).
Then, the selection is made by the following command:
dom.Particles[i]->Constitutive = 11;
Once the constitutive model for fluid is chosen, some parameters must be modified by a double value, the dynamic viscosity Mu for the Newtonian fluid. When the fluid is non-Newtonian, other parameters must be modified, such as reference dynamic viscosity MuRef and the yield stress T0, which are used in all the non-Newtonian fluid models. The regularisation parameter m used in Bingham and
Herschel-Bulkley models, which value is 300.0 by default; and time parameter t0mu used in Herschel-Bulkley. Thus, following the syntaxis can be implemented to modify the parameters depending on the selected model as shown below:
dom.Particles[i]->Mu = 1.0e-3;
dom.Particles[i]->MuRef = 2.0;
dom.Particles[i]->T0 = 5.0;
dom.Particles[i]->m = 500.0;
dom.Particles[i]->t0mu = 5.0e-4;
Such values depend on the experimental measurement. For further details regarding the models, review the papers cited above.
Since the system is assumed to be weakly compressible, an equation of state (EOS) must be used.
- 0 → Nearly incompressible fluids (Cole, 1965, p. 38-44).
- 1 → Artificial EOS (Morris et al., 2000).
- 2 → (Morris et al., 1997).
Thus, the EOS is selected as:
dom.Particles[i]->PresEq = 0;
Also, EOS 0 and 1 has a background pressure that is 0.0 by default, and can be modified using the following command:
dom.Particles[i]->P0 = P0;
The speed of the sound has an essential role since it influences the stability of the solution and the CPU time. Although a value is assigned in this example, some expressions can be found in the literature to compute the initial Cs, as shown in Morris et al. (1997). Typically, the used value is not the physical one. This is done to obtain results in a reasonable computational time, ensuring that it will not affect the quality of the solution.
dom.Particles[i]->Cs = Cs;
In order to diminish excessive oscillations caused by the shock fronts, an artificial viscous term is introduced. The filter can be employed by the following commands:
dom.Particles[i]->Alpha = 0.05;
dom.Particles[i]->Beta = 0.05;
Alpha and Beta are doubles whose values depend on the material and case. Their value can be between 0.0 and 1.0. To activate this artificial viscosity, either Alpha or Beta should be greater than zero. It has been found that 1.0 is an appropriate value for Alpha and Beta when simulating soils and close to 0.05 for fluid simulations.
SPH particles are prone to agglomerate in an unphysical way when particles have negative pressure coming from the equation of state. Monaghan (2000) introduce a term in the momentum equation to mitigate such a problem, which works as an artificial pressure to give repulsion between particles. It can be implemented in code using following lines:
dom.Particles[i]->TI = 0.3;
dom.Particles[i]->TIn = 4.0;
dom.Particles[i]->TIInitDist = dx;
By default, TIn is 4.0 in the code as recommended by Gray et al. (2001) for solid materials. To activate this option in the code, TI should be greater than zero. Generally, TI and TIn are recommended to be 0.3 and 4.0 for solids and fluids (Gray et al., 2001; Monaghan 2000), 0.5 and 2.55 for soils with cohesion (Bui et al., 2008).
When a particle is near the boundary or to the free surface, the kernel function suffers from a lack of particles. A manner to correct such a lack is using the Shepard filter so that the density of the particles can be reinitialised. This filter is applied every M time-steps (usually M is on the order of 30 time-steps) (Crespo 2008). The following command must be used to activate this filter in the code:
dom.Particles[a]->Shepard = true;
If the number of neighbouring particles is not enough (using the number density concept), the Shepard filter will be skipped for that particle.
In summary, the physical and numerical parameters for the FLUID are written below:
dom.Particles[i]->Material = 1;
dom.Particles[i]->Constitutive = 11;
dom.Particles[i]->Mu = Mu;
dom.Particles[i]->MuRef = Mu;
dom.Particles[i]->PresEq = 0;
dom.Particles[i]->P0 = P0;
dom.Particles[i]->Cs = CsF;
dom.Particles[i]->Alpha = 0.05;
dom.Particles[i]->Beta = 0.05;
dom.Particles[i]->TI = 0.3;
dom.Particles[i]->TIn = 4.0;
dom.Particles[i]->TIInitDist = dx;
dom.Particles[i]->Shepard = true;
It is necessary to define the physical and numerical parameters for the soil as well. The Material can be chosen by setting up an integer of value 3. There are two options for the Constitutive model for soils, 31->Drucker-Prager and 32->Hypoplastic.
dom.Particles[i]->Material = 3;
dom.Particles[i]->Constitutive = 31;
The failure criterion must be chosen when selecting an elastoplastic model. There is only one option for solids and two options for soils.
- 1 → Von-Mises failure criterion
- 2 → Drucker-Prager failure criterion - associated flow rule
- 3 → Drucker-Prager failure criterion - non-associated flow rule
If Drucker-Prager failure criterion with a non-associated flow rule is used, the required option is the third:
dom.Particles[i]->Fail = 3;
The assignation of some physical parameters previously defined are shown below. In any of the models for soils, elastoplastic or hypoplastic, a cohesion c can be set up as
follows:
dom.Particles[i]->c = c;
Other parameters that are proper from elastoplastic models are assign based on their previous definition as shown above, for instance: angle of shearing resistance phi, dilatancy angle psi, bulk modulus K, shear modulus G.
dom.Particles[i]->phi = (Phi*M_PI)/180.0;
dom.Particles[i]->psi = (Psi*M_PI)/180.0;
dom.Particles[i]->K = K;
dom.Particles[i]->G = G;
If the hypoplastic model is chosen, the option for such a Constitutive model is 32. Four parameters have to be configured and will describe the material behaviour, in addition to the cohesion that was previously shown. The physical meaning of those four parameters have not been found up to now, being an underlined characteristic in most of the hypoplastic models. For example, the following values were taken from Peng et al. (2015):
dom.Particles[i]->c1hp = -58.3;
dom.Particles[i]->c2hp = -178.5;
dom.Particles[i]->c3hp = -589.0;
dom.Particles[i]->c4hp = -182.7;
Since there is an interaction between fluid and soil, a seepage force can be defined based on the expected flow type in the porous media.
- 0 → Darcy's law
- 1 → Darcy's law & Kozeny-Carman coefficients
- 2 → Forchheimer equation & Ergun coefficients
- 2 → Forchheimer equation & Den Adel coefficients
dom.Particles[i]->SeepageType = 0;
Some parameters are needed for the fluid-soil interaction, where the intrinsic permeability k is the unique parameter when Darcy is chosen. If the other seepage equations are used, two more parameters are needed, porosity n0 and effective diameter d.
dom.Particles[i]->k = k;
dom.Particles[i]->n0 = ns;
dom.Particles[i]->d = 2.5e-4;
Also, the code has an option to have a variable porosity, as detailed by Korzani et al. (2018). Although variable porosity is not implemented in this case, if the modeller needs to activate this option, the following command can be used:
dom.Particles[i]->VarPorosity = true;
Water density should be given to soil particles to define the soil's buoyant unit weight in the code. However, the configured case is dynamic. so that the buoyancy is deactivated by defining the density as zero:
dom.Particles[i]->RhoF = 0.0;
In order to diminish excessive oscillations caused by the shock fronts, an artificial viscous term is introduced. The filter can be implemente via the following commands. However, there is a difference in the values for soils regarding the used in the fluid phase.
dom.Particles[i]->Alpha = 1.0;
dom.Particles[i]->Beta = 1.0;
SPH particles are prone to agglomerate in an unphysical way when particles have negative pressure coming from the equation of state. Monaghan (2000) introduce a term in the momentum equation to mitigate such a problem, which works as an artificial pressure to give repulsion between particles. It can be implemented in code using following lines:
dom.Particles[i]->TI = 0.5;
dom.Particles[i]->TIn = 2.55;
dom.Particles[i]->TIInitDist = dx;
By default, TIn is 4.0 in the code as recommended by Gray et al. (2001) for solid materials. To activate this option in the code, TI should be greater than zero. Generally, TI and TIn are recommended to be 0.5 and 2.55 for soils with cohesion (Bui et al, 2008).
The Shepard filter is also activated for soil particles, so that:
dom.Particles[i]->Shepard = true;
In summary, the physical and numerical parameters for the SOIL are written below:
dom.Particles[i]->Material = 3;
dom.Particles[i]->Constitutive = 31;
dom.Particles[i]->Fail = 3;
dom.Particles[i]->c = c;
dom.Particles[i]->phi = (Phi*M_PI)/180.0;
dom.Particles[i]->psi = (Psi*M_PI)/180.0;
dom.Particles[i]->K = K;
dom.Particles[i]->G = G;
dom.Particles[i]->SeepageType = 0;
dom.Particles[i]->k = k;
dom.Particles[i]->n = ns;
dom.Particles[i]->RhoF = 0.0;
dom.Particles[i]->Cs = CsS;
dom.Particles[i]->Alpha = 1.0;
dom.Particles[i]->Beta = 1.0;
dom.Particles[i]->TI = 0.5;
dom.Particles[i]->TIn = 2.55;
dom.Particles[i]->TIInitDist = dx;
dom.Particles[i]->Shepard = true;
Rectangle and cube are the primary shapes to discretise the domain by SPH, as pointed out before. If the modeller requires reshaping the discretised zones, it can be done by deleting some SPH particles from specific zones using basic functions. It can be done by changing the initial tag of the particles through the variable called dom.Particles[i]->ID = value. Then, the particles with that new ID can be eliminated using a command from the code called dom.DelParticles(ID) as will be shown below. Since the case is two-dimensional, the equation of a circle is used to empty a specific zone to insert a DEM sphere later:
if ( sqrt( pow(xb-Cx,2.0) + pow(yb-Cy,2.0) )<(Rc+0.5*dx) ) dom.Particles[i]->ID = 0;
An integer of value zero is the new ID of the unwanted particles. Out of the loop used to assign the physical and numerical parameters, the unwanted SPH particles can be deleted by the following command.
dom.DelParticles(0);
Domain shape is always a rectangular box. For obtaining a box size, two points are required, front-left-bottom and rear-right-top corners. For determining these points, maximum and minimum of x, y and z should be found using particles coordinates. This code finds these points coordinates automatically. However, if particles move out of this box, they will be deleted. Therefore, domain size should be extended in the case that deformation beyond the calculated box is required. Following syntaxes should be used for this extension:
dom.DomMax = Vec3_t( L+2.0*dx, H+0.2*H, 0.0 );
dom.DomMin = Vec3_t( -L-2.0*dx, -2.0*dx, 0.0 );
Declaring the DEM Particles
There are three kinds of DEM particles that can be used so far in the SPH-DEM coupling. Segments in 2D cases, Planes in 3D cases and Spheres that can be used in both 2D and 3D cases. To set up the other boundary conditions of a tank, three DEM segments will be implemented, left and right walls and the bottom. Also, an DEM sphere is inserted to simulate a particle with a diameter much larger than the effective diameter of the sand. The following lines are used to declare the three DEM segments and the sphere:
dom.AddSegment(/*Tag*/ -1,/*PositionB*/Vec3_t(-dx,-dx,0.0),/*PositionT*/Vec3_t(-dx,H+0.1*H,0.0),/*Thickness*/dx,/*Density*/rhoS); // Left wall
dom.AddSegment(/*Tag*/ -2,/*PositionB*/Vec3_t(L+dx,-dx,0.0),/*PositionT*/Vec3_t(L+dx,H+0.1*H,0.0),/*Thickness*/dx,/*Density*/rhoS); // Right wall
dom.AddSegment(/*Tag*/ -3,/*PositionL*/Vec3_t(-dx,-dx,0.0),/*PositionR*/Vec3_t(L+dx,-dx,0.0),/*Thickness*/dx,/*Density*/rhoS); // Bottom wall
dom.AddSphere(/*Tag*/ -4, /*Centre*/Vec3_t(Cx,Cy,Cz),/*Radius*/Rc,/*Density*/rhoB); // Sphere
The first parameter is a "Tag" to distinguish different subgroups of particles or even single particles. The vector Vec3_t(-dx,-dx,0.0) indicates the "Position L" of the left extreme of the DEM segment with coordinates Vec3_t(x, y, z). The vector Vec3_t(-dx,H+0.1*H,0.0) indicates the "Position R" of the right extreme of the DEM segment which coordinates are Vec3_t(x, y, z). "Thickness" is a double that defines the projection of the 2D segment into the third direction, z. However, it does not affect the simulation; a value equals dx is recommended. Then, one more parameter, "Density", is assigned; the same density used for the SPH particles is suggested. When more materials are declared into the domain, the maximum density of one of the components is suggested to set up the boundary conditions.
If the case is three-dimensional, planes are required. To add planes the following structure can be implemented:
dom.AddPlane(/*Tag*/ -5, /*Centre*/Vec3_t(Cx,Cy,Cz), /*Thickness*/dx, /*Length_x*/Lx, /*Length_y*/Ly, /*Density*/rhoS, /*Angle*/(-90.*M_PI)/180.0, /*Axis of rotation*/OrthoSys::e2 );
Vector Vec3_t(Cx,Cy,Cz) indicates the centre of the DEM particle, either sphere or plane. Lx and Ly indicate the length of the plane in each direction, x and y, respectively. The Angle is a double that specifies the angle through which the plane is rotated about an axis in radians. Axis of rotation indicates the axis about which the plane will be rotated. For example, OrthoSys::e0, OrthoSys::e1 or OrthoSys::e2 to rotate on x, y or z-axis, respectively.
After declaring the DEM particles, which are segments, in this case, it is necessary to assign some parameters to each DEM particle that define their role, to say, if they are boundary particles or free-moving objects. This will be done by using a loop for, which structure is written as:
for (size_t j = 0;j<dom.DEMParticles.Size();j++){
...
}
The gravity can be activated for DEM particles by setting up the force based on the DEM particle's mass and the gravity defined for the whole domain, as shown above.
dom.DEMParticles[j]->Ff = dom.DEMParticles[j]->Props.m*dom.Gravity;
Some of the properties of the DEM particles can be modified using Props. For example, if the SPH particle corresponds to the soil, then a value greater than zero should be used in the coefficient of friction dom.DEMParticles->Props.Mu of the DEM particles.
dom.DEMParticles[j]->Props.Mu = tan((Phi*M_PI)/180.0);
The lines shown below involve the computation of the minimum distance between an SPH particle and a DEM object based on the faces of the segment or plane, for instance. However, the sphere does not have such property, so that lines that involve Face of the DEM particles must be avoided. It can be done by the following command:
if ( dom.DEMParticles[j]->Tag == -4 ){continue;}
To set up the DEM particles as the boundary conditions, they have to be fixed. Hence, the following command is used for such a purpose:
dom.DEMParticles[j]->FixVeloc();
SPH particles do not have a physical radius, whereas DEM particles have. Besides, DEM particles have a halo to allow SPH particles to penetrate a certain distance into such a halo but no DEM particles per se (Trujillo-Vela et al., 2020). First, the thickness of the halo is defined as zero:
dom.DEMParticles[j]->Props.eps = 0.0;
Then, using a loop, the thickness of the halo of each DEM particle is modified based on the distance between the SPH and DEM particles as follows:
for (size_t i=0; i<dom.Particles.Size(); i++){
double dist = DEM::Distance(dom.Particles[i]->x,*dom.DEMParticles[j]->Faces[0]);
if (dist<dom.DEMParticles[j]->Props.eps || dom.DEMParticles[j]->Props.eps == 0.0) dom.DEMParticles[j]->Props.eps = dist;
}
Remark: When the DEM segments have an inclination and the SPH particles are placed at different distances from the DEM particles, the minimum distance can be used to define the thickness of the halo.
If some properties of specific DEM particles are needed to change, it can be performed as shown below.
dom.DEMParticles[0]->Props.Mu = 0.0;
dom.DEMParticles[3]->Props.eps = 0.5*dx;
Also, the stiffness of the DEM objects can be modified as follows:
dom.DEMstiff = 0.002;
Remark: DEMstiff = 0.1 by default. Nonetheless, in laboratory scale cases where the interaction forces between SPH particles and DEM objects are low, such value might produce some noise in the interaction. Thus, if the DEM stiffness is lower than 0.1, such noise can be reduced or avoided. Also, DEMstiff cannot be lower than 1.0e-4. Otherwise, SPH particles might go through the walls.
Starting the simulation
The following lines were introduced in order to check the value of some variables on the terminal:
std::cout<<"CsF = " <<CsF<<std::endl;
std::cout<<"CsS = " <<CsS<<std::endl;
std::cout<<"Steps = " <<steps<<std::endl;
std::cout<<"Resolution Sphere-Fluid = "<<(2.0*Rc/dx)<<std::endl;
The simulation can be started with the command Solve, so that:
dom.Solve(/*tf*/Tfinal,/*dt*/timestep,/*dtOut*/dtoutput,file_name,numfiles);
return 0;
The system will evolve in this case up to Tfinal = 2.0 seconds of simulation time. The next argument is the integration time step, timestep. The smaller, the more accurate the simulation is, and the more time-consuming. The 3rd argument is the report time-step, in other words, the time span between movie frames, previously defined as dtoutput. The next argument is the name of the output files, file_name. Finally, the last argument is the number of total files numfiles. These variables were defined above.
The way to execute the compiled file on the terminal is ./test_debris 4 &. Where 4 indicates the number of processors when required to change it from the maximum as it was explained before.
Simulation visualisation
MechSys uses VisIt to visualise simulation results. A full tutorial can be found here. After following the instructions and choosing a Pseudocolor->SPHCenter->Tag data type to be visualised, the dynamics of the system should resemble the following video:
Publications
- Korzani, M. G., Galindo Torres, S., Scheuermann, A., & Williams, D. J. (2016). Smoothed particle hydrodynamics into the fluid dynamics of classical problems. In Applied Mechanics and Materials (Vol. 846, pp. 73-78). Trans Tech Publications Ltd.
- Korzani, M. G., Galindo-Torres, S. A., Scheuermann, A., & Williams, D. J. (2017). Parametric study on smoothed particle hydrodynamics for accurate determination of drag coefficient for a circular cylinder. Water Science and Engineering, 10(2), 143-153.
- Korzani, M. G., Galindo-Torres, S. A., Scheuermann, A., & Williams, D. J. (2018). Smoothed Particle Hydrodynamics for investigating hydraulic and mechanical behaviour of an embankment under action of flooding and overburden loads. Computers and Geotechnics, 94, 31-45.
- Korzani, M. G., Galindo-Torres, S. A., Scheuermann, A., & Williams, D. J. (2018). SPH approach for simulating hydro-mechanical processes with large deformations and variable permeabilities. Acta Geotechnica, 13(2), 303-316.
- Trujillo-Vela, M. G., Galindo-Torres, S. A., Zhang, X., Ramos-Cañón, A. M., & Escobar-Vargas, J. A. (2020). Smooth particle hydrodynamics and discrete element method coupling scheme for the simulation of debris flows. Computers and Geotechnics, 125, 103669.
References
- Braune, L. & Lewiner, T. (2013). An initiation to sph. Rio de Janeiro: PUC-Rio, pages 1–7.
- Bui, H. H., Sako, K., & Fukagawa, R. (2007). Numerical simulation of soil–water interaction using smoothed particle hydrodynamics (SPH) method. Journal of Terramechanics, 44(5), 339-346.
- Bui, H. H., Fukagawa, R., Sako, K., & Ohno, S. (2008). Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic–plastic soil constitutive model. International journal for numerical and analytical methods in geomechanics, 32(12), 1537-1570.
- Bui, H. H., & Nguyen, G. D. Smoothed particle hydrodynamics (SPH) and its applications in geomechanics. ALERT Doctoral School 2020 Point based numerical methods in geomechanics, 3.
- Cole, R. H., & Weller, R. (1948). Underwater explosions. PhT, 1(6), 35.
- Dehnen, W., & Aly, H. (2012). Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Notices of the Royal Astronomical Society, 425(2), 1068-1082.
- Ghadimi, P., Farsi, M., and Dashtimanesh, A. (2012). Study of various numerical aspects of 3d-sph for simulation of the dam break problem. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 34(4):486–491.
- Gray, J. P., Monaghan, J. J., & Swift, R. P. (2001). SPH elastic dynamics. Computer methods in applied mechanics and engineering, 190(49-50), 6641-6662.
- Ké, D. D., & Turcotte, G. (1980). Viscosity of biomaterials. Chemical Engineering Communications, 6(4-5), 273-282.
- Monaghan, J. J. (2000). Sph without a tensile instability. Journal of Computational Physics, 159(2):290–311.
- Monaghan, J. J. (2005). Smoothed particle hydrodynamics. Reports on progress in physics, 68(8):1703–1759.
- Morris, J. P., Fox, P. J., & Zhu, Y. (1997). Modeling low Reynolds number incompressible flows using SPH. Journal of computational physics, 136(1), 214-226.
- Morris, J. P. (2000). Simulating surface tension with smoothed particle hydrodynamics. International journal for numerical methods in fluids, 33(3), 333-353.
- Nguyen, C. T., Nguyen, C. T., Bui, H. H., Nguyen, G. D., & Fukagawa, R. (2017). A new SPH-based approach to simulation of granular flows using viscous damping and stress regularisation. Landslides, 14(1), 69-81.
- Imran, J., Parker, G., Locat, J., & Lee, H. (2001). 1D numerical model of muddy subaqueous and subaerial debris flows. Journal of hydraulic engineering, 127(11), 959-968.
- Papanastasiou, T. C. (1987). Flows of materials with yield. Journal of Rheology, 31(5), 385-404.
- Peng, C., Wu, W., Yu, H. S., & Wang, C. (2015). A SPH approach for large deformation analysis with hypoplastic constitutive model. Acta Geotechnica, 10(6), 703-717.
- Shao, S., & Lo, E. Y. (2003). Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Advances in water resources, 26(7), 787-800.
- Takeda, H., Miyama, S. M., & Sekiya, M. (1994). Numerical simulation of viscous flow by smoothed particle hydrodynamics. Progress of theoretical physics, 92(5), 939-960.