Smoothed particle hydrodynamics

Introduction

The example shown in Figure 1 contrasts the two approaches. Both discrete representations model the initial configuration of a body of fluid inside a bottle, as described in detail in Impact of a water-filled bottle. The model on the left is a traditional tetrahedron mesh of the volume occupied by the fluid. On the right, the same volume is represented by a collection of discrete points. Note that in the latter case there are no edges connecting the points as these points (pseudo-particles) do not require the definition of multiple-node element connectivity, as is the case in the traditional finite element representation on the left. An alternative to directly defining particle elements is to define conventional continuum finite elements and automatically convert them to particle elements at the beginning of the analysis or during the analysis, as discussed in Finite element conversion to SPH particles.

Figure 1. Finite element mesh and SPH particle distribution.

Smoothed particle hydrodynamics is a fully Lagrangian modeling scheme permitting the discretization of a prescribed set of continuum equations by interpolating the properties directly at a discrete set of points distributed over the solution domain without the need to define a spatial mesh. The method's Lagrangian nature, associated with the absence of a fixed mesh, is its main strength. Difficulties associated with fluid flow and structural problems involving large deformations and free surfaces are resolved in a relatively natural way.

At its core, the method is not based on discrete particles (spheres) colliding with each other in compression or exhibiting cohesive-like behavior in tension as the word particle might suggest. Rather, it is simply a clever discretization method of continuum partial differential equations. In that respect, smoothed particle hydrodynamics is quite similar to the finite element method. SPH uses an evolving interpolation scheme to approximate a field variable at any point in a domain. The value of a variable at a particle of interest can be approximated by summing the contributions from a set of neighboring particles, denoted by subscript j, for which the “kernel” function, W, is not zero

⟨ f (x) ⟩ ≃ \sum_{j} \frac{m_{j}}{ρ_{j}} f_{j} W (| x - x_{j} |, h) .

An example kernel function is shown in Figure 2. The smoothing length, h, determines how many particles influence the interpolation for a particular point.

Figure 2. Kernel function.

The SPH method has received substantial theoretical support since its inception (Gingold and Monaghan, 1977), and the number of publications related to the method is now very large. A number of references are listed below.

The method can use any of the materials available in Abaqus/Explicit (including user materials). You can specify initial conditions and boundary conditions as for any other Lagrangian model. Contact interactions with other Lagrangian bodies are also allowed, thus expanding the range of applications for which this method can be used.

The method is less accurate in general than Lagrangian finite element analyses when the deformation is not too severe and than coupled Eulerian-Lagrangian analyses in higher deformation regimes. If a large percentage of all nodes in the model are associated with smoothed particle hydrodynamics, the analysis may not scale well if multiple CPUs are used.

Applications

Smoothed particle hydrodynamic analyses are effective for applications involving extreme deformation. Fluid sloshing, wave engineering, ballistics, spraying (as in paint spraying), gas flow, and obliteration and fragmentation followed by secondary impacts are a few examples. There are many applications for which both the coupled Eulerian-Lagrangian and the smoothed particle hydrodynamic methods can be used. In many coupled Eulerian-Lagrangian analyses the material to void ratio is small and, consequently, the computational effort may be prohibitively high. In these cases, the smoothed particle hydrodynamic method is preferred. For example, tracking fragments from primary impacts through a large volume until secondary impact occurs can be very expensive in a coupled Eulerian-Lagrangian analysis but comes at no additional cost in a smoothed particle hydrodynamic analysis.

Impact of a water-filled bottle includes an example of using the smoothed particle hydrodynamic method to model the violent sloshing associated with the impact.

Artificial viscosity

Artificial viscosity in smoothed particle hydrodynamics has the same meaning as bulk viscosity for finite elements. Similar to other Lagrangian elements, particle elements use linear and quadratic viscous contributions to dampen high frequency noise from the computed response. In rare cases when the default values are not appropriate, you can control the amount of artificial viscosity included in a smoothed particle hydrodynamic analysis.

Input File Usage

Use the following options to specify scale factors for the linear and quadratic artificial viscosities:

SECTION CONTROLS
 , , , scale factor for linear artificial viscosity, scale factor for quadratic artificial viscosity

SPH tensile instability control

The classical SPH formulation can display a tensile instability that results in an unphysical clustering of particles. The clustering is particularly noticeable when materials are in a tensile stress state. To relieve the SPH tensile instability, Abaqus uses the artificial stress approach of Monaghan, 2000. In this method an artificial stress contribution that acts as a short-range repulsion force is added to mitigate the tensile instability. This repulsion force is used to prevent particles from getting too close to each other. This force is related to the inter-particle distance: the smaller the distance the greater the force. This feature is effective only when an initially uniform distribution of particles is used.

Input File Usage

Use the following option to specify the SPH tensile instability control:

SECTION CONTROLS, SPH TENSILE INSTABILITY CONTROL=YES

Initial conditions

Initial conditions in Abaqus/Standard and Abaqus/Explicit describes all of the initial conditions that are available for an explicit dynamic analysis. Initial conditions pertinent to mechanical analyses can be used in a smoothed particle hydrodynamic analysis.

Boundary conditions

Boundary conditions are defined as described in Boundary conditions in Abaqus/Standard and Abaqus/Explicit. The default interaction between SPH particles and Lagrangian surfaces is based on node-to-surface contact. To improve the SPH solution near the Lagrangian boundaries for flow-type problems, you can define the surface behavior explicitly.

Improving the SPH boundary behavior near Lagrangian boundaries

In the SPH method, particles interact with their neighboring particles through a smoothing kernel function. The smoothing kernel function has a finite supporting domain. To achieve accurate SPH simulation results, there should be enough neighboring particles (50–60) within the supporting domain of the smoothing kernel. However, special considerations are needed near the boundaries. When particles are close to Lagrangian boundaries, part of the supporting domain of the smoothing kernel (that beyond the Lagrangian boundaries) will not be filled with any SPH particles. Therefore, the integration accuracy for those SPH particles close to the Lagrangian boundaries may be negatively affected. In addition, the disturbance of the solution close to the boundaries may eventually lead to inaccurate simulation results for the whole field.

To improve the SPH solution near the Lagrangian boundaries, you can specify surface behavior that is based on the ghost particle method (Colagrossi and Landrini, 1977) and define the tangential-type of SPH boundary behavior. When SPH particles are close to Lagrangian boundary surfaces, a virtual plane is formed by taking the average of the adjacent boundary surfaces; based on this virtual plane, ghost particles are created to help improve integration accuracy. The state (pressure, density, velocity, etc.) of the ghost particles is computed from the real physical particles. With the improved boundary condition, SPH particles on opposite sides of a surface do not interact with each other.

Input File Usage

Use one of the following options to specify the SPH surface behavior:

SPH SURFACE BEHAVIOR, TYPE=FREE SLIP (default) 
SPH SURFACE BEHAVIOR, TYPE=NO SLIP

Loads

The loading types available for an explicit dynamic analysis are explained in About loads. Concentrated nodal loads can be applied as usual. Gravity loads are the only distributed loads allowed in a smoothed particle hydrodynamic analysis.

Smoothed particle hydrodynamic domain

The region within which the particle search (finding all neighbors for all particles) is performed is discussed in Using section controls to define the particle tracking box for DEM and SPH particles.

Material options

Any of the material models in Abaqus/Explicit can be used in a smoothed particle hydrodynamic analysis.

Elements

The smoothed particle hydrodynamic method is implemented via the formulation associated with PC3D elements. These 1-node elements are simply a means of defining particles in space that model a particular body or bodies. These particle elements utilize existing functionality in Abaqus to reference element-related features such as materials, initial conditions, distributed loads, and visualization.

You define these elements in a similar fashion as you would define point masses. The coordinates of these points lie either on the surface or in the interior of the body being modeled, similar to the nodes of a body meshed with brick elements. For more accurate results, you should strive to space the nodal coordinates of these particles as uniformly as possible in all directions.

An alternative to directly defining PC3D elements is to define conventional continuum finite element types C3D8R, C3D6, or C3D4 and automatically convert them to particle elements at the beginning of the analysis or during the analysis, as discussed in Finite element conversion to SPH particles.

The smoothed particle hydrodynamic method implemented in Abaqus/Explicit uses a cubic spline as the interpolation polynomial and is based on the classical smoothed particle hydrodynamic theory as outlined in the references below.

The smoothed particle hydrodynamic method is not implemented for two-dimensional elements. Axisymmetry can be simulated using a wedge-shaped sector and symmetry boundary conditions. There are no hourglass or distortion control forces associated with PC3D elements. These elements do not have faces or edges associated with them.

SPH kernel interpolator and formulations

By default, the smoothed particle hydrodynamic method implemented in Abaqus/Explicit uses a cubic spline as the interpolation polynomial. Alternatively, you can choose a quadratic (Johnson et al, 1996) or quintic (Wendland, 1995) interpolator.

The implementation is based on the classical smoothed particle hydrodynamic theory as outlined in the references below. You also have the option of using a mean flow correction configuration update, commonly referred to in the literature as the XSPH method (see Monaghan, 1992), as well as the corrected kernel of Randles and Libersky, 1997, also referred to as the normalized SPH (NSPH) method. In the XSPH method, the mean velocity filtering coefficient is required for modifying the coordinate update for the particles. A zero value for this coefficient leads to the classical SPH method.

You can control these settings as discussed in Using section controls for smoothed particle hydrodynamics (SPH).

Computing the particle volume

There is currently no capability to automatically compute the volume associated with these particles. Hence, you need to supply a characteristic length that will be used to compute the particle volume, which in turn is used to compute the mass associated with the particle. It is assumed that the nodes are distributed uniformly in space and that each particle is associated with a small cube centered at the particle. When stacked together, these cubes will fill the overall volume of the body with some minor approximation at the free surface of the body. The characteristic length is half the length of the cube side. From a practical perspective, once you have created the nodes, you can use the half-distance between two nodes as the characteristic length. Alternatively, if you know the mass and density of the part, you can compute the volume of the part and divide it by the total number of particles in the part to obtain the volume of the small cube associated with each particle. Half of the cubic root of this small volume is a reasonable characteristic length for this particle set. You can check the mass of individual sets in the model if you request that model definition data be printed to the data (.dat) file (see Model and history definition summaries).

Input File Usage

Use the following options to define a smoothed particle hydrodynamic body:

ELEMENT, TYPE=PC3D, ELSET=particle_body
element number, node number

Repeat the data line as often as necessary.

SOLID SECTION, ELSET=particle_body, MATERIAL=material_name
characteristic length associated with particle volume

Smoothing length calculation

Even though particle elements are defined in the model using one node per element, the smoothed particle hydrodynamic method computes contributions for each element based on neighboring particles that are within a sphere of influence. The radius of this sphere of influence is referred to in the literature as the smoothing length. The smoothing length is independent of the characteristic length discussed above and governs the interpolation properties of the method. By default, the smoothing length is computed automatically. As the deformation progresses, particles move with respect to each other and, hence, the neighbors of a given particle can (and typically do) change. Every increment Abaqus/Explicit recomputes this local connectivity internally and computes kinematic quantities (such as normal and shear strains, deformation gradients, etc.) based on contributions from this cloud of particles centered at the particle of interest. Stresses are then computed in a similar fashion as for reduced-integration brick elements, which are in turn used to compute element nodal forces for the particles in the cloud based on the smoothed particle hydrodynamic formulation.

By default, Abaqus/Explicit computes a smoothing length at the beginning of the analysis such that the average number of particles associated with an element is roughly between 50 and 60. The smoothing length can be kept constant or can vary during the analysis. In both cases the average number of particles per element can either decrease or increase during the analysis depending on whether the average behavior in the model is expansive or compressive, respectively.

If the average number of particles per element decreases below a certain minimum number, the deformation gradient may not be calculated accurately due to the limited support from surrounding particles within the sphere of influence. In this case the deformation gradient is kept constant between the previous time increment and the current time increment. By default, the minimum number of particles associated with one element is 5.

You can control most of these settings as discussed in Using section controls for smoothed particle hydrodynamics (SPH).

Smoothed particle hydrodynamic domain

By default, all particles are tracked as the analysis progresses. Alternatively, you can define a particle tracking box to define the region within which the particle search (finding all neighbors for all particles) is performed as discussed in Using section controls to define the particle tracking box for DEM and SPH particles.

Constraints

Since the PC3D elements are Lagrangian elements, their nodes can be involved in other features, such as other elements, connectors, or constraints. Since these elements do not have faces or edges, an element-based surface cannot be defined using PC3D elements. Consequently, constraints that require element-based surfaces (such as fasteners) cannot be defined for particles.

Interactions

Bodies modeled with particles can interact with other finite element meshed bodies via contact. The contact interaction is the same as any contact interaction between a node-based surface (associated with the particles) and an element-based or analytical surface. Both general contact and contact pairs can be used. All interaction types and formulations available for contact involving a node-based surface are allowed, including cohesive behavior. Different contact properties can be assigned via the usual options. By default, the particles are not part of the general contact domain similar to other 1-node elements (such as point masses). The default contact thickness for particles is the same value specified as the characteristic length on the section definition, and it is constant through the analysis. Thus, for contact purposes, particles behave as spheres with radii equal to the radius of a sphere inscribed in the small cube associated with the particle volume as described above.

You should not specify a contact thickness of zero for the nodes associated with PC3D elements or contact may not be resolved robustly. The recommended approach is to use the default or specify a reasonable contact thickness.

Interaction between different bodies all modeled with PC3D elements is allowed. However, this interaction is meaningful only in cases when the colliding smoothed particle hydrodynamic bodies are made of the same fluid-like material, such as a water drop falling in a bucket partially filled with water. In solids-related applications, such as modeling a bullet perforating an armor plate, one of the bodies must be modeled using regular finite elements.

Contact interactions cannot be defined between particles and Eulerian regions.

Input File Usage

Use the following options to define contact between a meshed or analytical surface with a particle-based surface:

CONTACT
CONTACT INCLUSIONS
node-based particle surface, element-based/analytical_surface

Output

The element output available for PC3D elements includes all mechanics-related output for continuum elements: stress; strain; energies; and the values of state, field, and user-defined variables. The nodal output includes all output variables generally available in Abaqus/Explicit analyses.

Particles can be visualized in Abaqus/CAE via circular discs. In contour plots the values of field output variables are shown as circular patches of color. Symbol plots are also available.

Limitations

Smoothed particle hydrodynamic analyses are subject to the following limitations:

They are less accurate in general than Lagrangian finite element analyses when the deformation is not too severe and the elements are not distorted. In higher deformation regimes coupled Eulerian-Lagrangian analyses are also generally more accurate. The smoothed particle dynamic method should be used primarily in cases when the conventional finite element method or the coupled Eulerian-Lagrangian method have reached their inherent limitations or are prohibitively expensive to perform.
When the material is in a state of tensile stress, the particle motion may become unstable leading to the so-called tensile instability. This instability, which is strictly related to the interpolation technique of the standard smoothed particle dynamic method, is especially noticeable when simulating the stretched state of a solid. As a consequence, particles tend to clump together and show fracture-like behavior.
Mass distribution in a body defined with particle elements is somewhat different when compared to the mass distribution of the same body defined with continuum elements, such as C3D8R elements. When particle elements are used, the volume of all particles in that body are the same. Consequently, the nodal mass associated with all particles in that body is the same. If the nodes are not placed in a regular cubic arrangement, the mass distribution is somewhat inexact, particularly at the free surface of the body being modeled.
Surface loads cannot be specified on PC3D elements. However, distributed loads, such as pressure, can be applied to other finite element surfaces that can apply a pressure onto the particle elements via contact interactions.
Bodies modeled with particles that were not defined using the same section definition do not interact with each other. Hence, you cannot use smoothed particle hydrodynamics to model the mixing of bodies with dissimilar materials.
Selective subcycling and the contact pair algorithm are not supported.
Rough friction should not be used on surface interactions defined for SPH particles.
The functionality is not directly supported in Abaqus/CAE. However, you can do the following:
- You can use the existing functionality in Abaqus/CAE to generate mass elements, write an input file, and then manually edit the input file to convert the mass elements to particles.
- You can use finite element conversion to SPH particles (see Boundary conditions) at the beginning of an analysis (a time-based conversion criterion at time equal to zero).
- You can create a mesh using C3D8R elements, write an input file, and then use a script to convert these elements to SPH particles. An example of a Python script to convert solid elements to SPH particles is described in “Generating SPH particle elements from a solid mesh” in the Dassault Systèmes Knowledge Base at http://www.3ds.com/support/knowledge-base.
Within a given body (part) defined via one solid section definition, gravity loads and mass scaling cannot be specified selectively on a subset of elements referenced by this definition. Instead, the two features must be applied to all the elements in the element set associated with the solid section definition.

Smoothed particle hydrodynamic computations are distributed across parallel domains in most cases; however, they are all performed by a single domain (with a single processor) for models with any of the following characteristics (which will often dramatically degrade parallel scalability):

Finite element conversion to SPH particles (see Boundary conditions) after the beginning of an analysis. The parallel smoothed particle hydrodynamics implementation supports conversion to SPH particles only at the beginning of an analysis (at time equal to zero).
Multiple solid sections for PC3D elements
Normalized kernel specified as a section control
Predefined field variable (including temperature) dependence of material properties

Smoothed particle hydrodynamic analyses are subject to the following limitations if multiple CPUs are used:

Contact output is not supported for smoothed particle hydrodynamic slave nodes.
Element history output is not supported.
Dynamic load balancing cannot be activated.
If any SPH particles participate in general contact, all SPH particles must be included in the general contact definition.
At least 10,000 particles per domain is suggested to achieve good scalability.
A significant increase in memory usage may be needed if a large number of CPUs are used.

Input file template

The following example illustrates a smoothed particle hydrodynamic analysis of a bottle filled with fluid being dropped on the floor. The plastic bottle and the floor are modeled with conventional shell elements. The fluid is modeled via smoothed particle hydrodynamics using PC3D elements. The nodal coordinates of the particles are defined such they are all located inside the bottle. Material property definitions are defined as usual for both the fluid and the bottle. Contact interaction is defined between the smoothed particle hydrodynamic particles representing the water (node-based surface) and the interior walls of the bottle and also between the bottle exterior and the floor using element-based surfaces (not shown). Output is requested for stresses (pressure) and density in the fluid.

HEADING
…
ELEMENT, TYPE=PC3D, ELSET=Fluid_Inside_The_Bottle
Element number, node number
…
SOLID SECTION, ELSET=Fluid_Inside_The_Bottle, MATERIAL=Water
Element characteristic length associated with particle volume
MATERIAL, NAME=Water
Material definition for water, such as an EOS material
ELEMENT, TYPE=S4R, ELSET=Plastic_Bottle
Element definitions for the shells
**
INITIAL CONDITIONS, TYPE=VELOCITY
Data lines to define velocity initial conditions
NSET, NSET=Water_Nodes, ELSET=Fluid_Inside_The_Bottle
SURFACE, NAME=Water_Surface, TYPE=NODE
Water_Nodes,
SURFACE, NAME=Bottle_Interior
Plastic_Bottle, SNEG
**
STEP
DYNAMIC, EXPLICIT
DLOAD
Data lines to define gravity load
**
CONTACT
CONTACT INCLUSIONS
Water_Surface, Bottle_Interior
**
OUTPUT, FIELD
ELEMENT OUTPUT, ELSET=Fluid_Inside_The_Bottle
S, DENSITY
END STEP

References

Colagrossi, A., and M. Landrini, “Numerical Simulation of Interfacial Flows by Smoothed Particle Hydrodynamics,” Journal of Computational Physics, 2003.
Gingold, R. A., and J. J. Monaghan, “Smoothed Particle Hydrodynamics: Theory and Application to Non-Spherical Stars,” Royal Astronomical Society, Monthly Notices, vol. 181, pp. 375–389, 1977.
Johnson, J., R. Stryk, and S. Beissel, “SPH for High Velocity Impact Calculations,” Computer Methods in Applied Mechanics and Engineering, 1996.
Libersky, L. D., and A. G. Petschek, “High Strain Lagrangian Hydrodynamics,” Journal of Computational Physics, vol. 109, pp. 67–75, 1993.
Monaghan, J., “Smoothed Particle Hydrodynamics,” Annual Review of Astronomy and Astrophysics, 1992.
Monaghan, J., “SPH without a Tensile Instability,” Journal of Computational Physics, vol. 159, 2000.
Munjiza, A., and K. R. F. Andrews, “NBS Contact Detection Algorithm for Bodies of Similar Size,” International Journal for Numerical Methods in Engineering, vol. 43, pp. 131–149, 1998.
Randles, P. W., and L. D. Libersky, “Recent Improvements in SPH Modeling of Hypervelocity Impact,” International Journal of Impact Engineering, 1997.
Swegle, J. W., and S. W. Attaway, “An Analysis of Smoothed Particle Hydrodynamics,” Sandia National Lab Report SAND93–2513, 1994.
Wendland, H., “Piecewise Polynomial, Positive Definite and Compactly Supported Radial Functions of Minimal Degree,” Advances in Computational Mathematics, 1995.