Bases: object
Bases: pysph.sph.equation.Equation
Add a body force to the particles:
\(\boldsymbol{f} = f_x, f_y, f_z\)
Bases: pysph.sph.equation.Equation
Density rate:
\(\frac{d\rho_a}{dt} = \sum_b m_b \boldsymbol{v}_{ab}\cdot \nabla_a W_{ab}\)
Bases: pysph.sph.equation.Equation
Compute the pressure using the Isothermal equation of state:
\(p = p_0 + c_0^2(\rho_0 - \rho)\)
Bases: pysph.sph.equation.Equation
Classical Monaghan style artificial viscosity [Monaghan2005]
where
with
[Monaghan2005] | J. Monaghan, “Smoothed particle hydrodynamics”, Reports on Progress in Physics, 68 (2005), pp. 1703-1759. |
Bases: pysph.sph.equation.Equation
Good old Summation density:
\(\rho_a = \sum_b m_b W_{ab}\)
Bases: pysph.sph.equation.Equation
Compute the SPH evaluation for the velocity gradient tensor in 2D.
The expression for the velocity gradient is:
\(\frac{\partial v^i}{\partial x^j} = \sum_{b}\frac{m_b}{\rho_b}(v_b - v_a)\frac{\partial W_{ab}}{\partial x_a^j}\)
The tensor properties are stored in the variables v_ij where ‘i’ refers to the velocity component and ‘j’ refers to the spatial component. Thus v_21 is \(\frac{\partial v}{\partial x}\)
Bases: pysph.sph.equation.Equation
Position stepping with XSPH correction [Monaghan1992]
[Monaghan1992] | J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574. |
This equation must be used to advect the particles. XSPH can be turned off by setting the parameter eps = 0.
Bases: pysph.sph.equation.Equation
The XSPH correction [Monaghan1992] alone. This is meant to be used with a leap-frog integrator which already considers the velocity of the particles. It simply computes the correction term and adds that to ax, ay, az.
[Monaghan1992] | J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574. |
This equation must be used to advect the particles. XSPH can be turned off by setting the parameter eps = 0.
Bases: pysph.sph.equation.Equation
Continuity equation with dissipative terms
\(\frac{d\rho_a}{dt} = \sum_b \rho_a \frac{m_b}{\rho_b} \left( \boldsymbol{v}_{ab}\cdot \nabla_a W_{ab} + \delta \eta_{ab} \cdot \nabla_{a} W_{ab} (h_{ab}\frac{c_{ab}}{\rho_a}(\rho_b - \rho_a)) \right)\)
[Marrone2011] | S. Marrone et al., “delta-SPH model for simulating violent impact flows”, Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp 1526–1542. |
Bases: pysph.sph.equation.Equation
Classic Monaghan Style Momentum Equation with Artificial Viscosity
where
with
[Monaghan1992] | J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574. |
Bases: pysph.sph.equation.Equation
Momentum equation defined in JOSEPHINE and the delta-SPH model
where
[Marrone2011] | S. Marrone et al., “delta-SPH model for simulating violent impact flows”, Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp 1526–1542. |
[Cherfils2012] | J. M. Cherfils et al., “JOSEPHINE: A parallel SPH code for free-surface flows”, Computer Physics Communications, 183 (2012), pp 1468–1480. |
Artificial viscosity is used in this momentum equation and is controlled by the parameter \(\alpha\). This form of the artificial viscosity is similar but not identical to the Monaghan-style artificial viscosity.
Bases: pysph.sph.equation.Equation
Pressure gradient discretized using number density:
Bases: pysph.sph.equation.Equation
Tait equation of state for water-like fluids
\(p_a = \frac{c_{0}^2\rho_0}{\gamma}\left( \left(\frac{\rho_a}{\rho_0}\right)^{\gamma} -1\right)\)
[Cole1948] | H. R. Cole, “Underwater Explosions”, Princeton University Press, 1948. |
[Batchelor2002] | G. Batchelor, “An Introduction to Fluid Dynamics”, Cambridge University Press, 2002. |
[Monaghan2005] | J. Monaghan, “Smoothed particle hydrodynamics”, Reports on Progress in Physics, 68 (2005), pp. 1703-1759. |
The reference speed of sound, c0, is to be taken approximately as 10 times the maximum expected velocity in the system. The particle sound speed is given by the usual expression:
\(c_a = \sqrt{\frac{\partial p}{\partial \rho}}\)
Bases: pysph.sph.equation.Equation
Tait Equation of State with Hughes and Graham Correction
where
[Hughes2010] | J. P. Hughes and D. I. Graham, “Comparison of incompressible and weakly-compressible SPH models for free-surface water flows”, Journal of Hydraulic Research, 48 (2010), pp. 105-117. |
The correction is to be applied on boundary particles and imposes a minimum value of the density (rho0) which is set upon instantiation. This correction avoids particle sticking behaviour at walls.
Bases: pysph.sph.equation.Equation
Update the particle smoothing lengths
\(h_a = hdx \left(\frac{m_a}{\rho_a}\right)^{\frac{1}{d}}\)
[Ferrari2009] | A. Ferrari et al., “A new 3D parallel SPH scheme for free surface flows”, Computers and Fluids, 38 (2009), pp. 1203–1217. |
Ideally, the kernel scaling factor should be determined from the kernel used based on a linear stability analysis. The default value of (hdx=1) reduces to the formulation suggested by Ferrari et al. who used a Cubic Spline kernel.
Typically, a change in the smoothing length should mean the neighbors are re-computed which in PySPH means the NNPS must be updated. This equation should therefore be placed as the last equation so that after the final corrector stage, the smoothing lengths are updated and the new NNPS data structure is computed.
Note however that since this is to be used with incompressible flow equations, the density variations are small and hence the smoothing lengths should also not vary too much.
Viscosity functions
Bases: pysph.sph.equation.Equation
Artificial viscosity proposed By P. Cleary:
\[\mathcal{Pi}_{ab} = -\]
rac{16}{mu_a mu_b}{ ho_a ho_b
(mu_a + mu_b)}left(
ight),
where the viscosity is determined from the parameter \(lpha\) as
\[\mu_a = \]
rac{1}{8}lpha h_a c_a ho_a
This equation is described in the 2005 review paper by Monaghan
- J. J. Monaghan, “Smoothed Particle Hydrodynamics”, Reports on Progress in Physics, 2005, 68, pp 1703–1759 [JM05]
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
[Adami2012] (1, 2) S. Adami et. al “A generalized wall boundary condition for smoothed particle hydrodynamics”, Journal of Computational Physics (2012), pp. 7057–7075.
[Adami2013] S. Adami et. al “A transport-velocity formulation for smoothed particle hydrodynamics”, Journal of Computational Physics (2013), pp. 292–307.
Bases: pysph.sph.equation.Equation
Conservation of mass equation
Eq (6) in [Adami2012]:
Bases: pysph.sph.equation.Equation
Artificial stress contribution to the Momentum Equation
where the artificial stress terms are given by:
Bases: pysph.sph.equation.Equation
Artificial viscosity for the momentum equation
Eq. (11) in [Adami2012]:
where
Bases: pysph.sph.equation.Equation
Momentum equation for the Transport Velocity Formulation: Pressure
Eq. (8) in [Adami2013]:
where
This equation should have the destination as fluid and sources as fluid and boundary particles.
This function also computes the contribution to the background pressure and accelerations due to a body force or gravity.
The body forces are damped according to Eq. (13) in [Adami2012] to avoid instantaneous accelerations. By default, damping is neglected.
Bases: pysph.sph.equation.Equation
Momentum equation for the Transport Velocity Formulation: Viscosity
Eq. (8) in [Adami2013]:
where
Bases: pysph.sph.equation.Equation
Extrapolating the fluid velocity on to the wall
Eq. (22) in [Adami2012]:
The destination particle array for this equation should define the filtered velocity variables \(uf, vf, wf\).
Bases: pysph.sph.equation.Equation
Solid wall boundary condition [Adami2012]
This boundary condition is to be used with fixed ghost particles in SPH simulations and is formulated for the general case of moving boundaries.
The velocity and pressure of the fluid particles is extrapolated to the ghost particles and these values are used in the equations of motion.
No-penetration:
Ghost particles participate in the continuity and state equations with fluid particles. This means as fluid particles approach the wall, the pressure of the ghost particles increases to generate a repulsion force that prevents particle penetration.
No-slip:
Extrapolation is used to set the dummy velocity of the ghost particles for viscous interaction. First, the smoothed velocity field of the fluid phase is extrapolated to the wall particles:
In the second step, for the viscous interaction in Eqs. (10) in [Adami2012] and Eq. (8) in [Adami2013], the velocity of the ghost particles is assigned as:
where \(v_w\) is the prescribed wall velocity and \(v_b\) is the ghost particle in the interaction.
For this equation the destination particle array should be the fluid and the source should be ghost or boundary particles. The boundary particles must define a prescribed velocity \(u_0, v_0, w_0\)
Bases: pysph.sph.equation.Equation
Solid wall pressure boundary condition [Adami2012]
This boundary condition is to be used with fixed ghost particles in SPH simulations and is formulated for the general case of moving boundaries.
The velocity and pressure of the fluid particles is extrapolated to the ghost particles and these values are used in the equations of motion.
Pressure boundary condition:
The pressure of the ghost particle is also calculated from the fluid particle by interpolation using:
where the subscripts g and f relate to the ghost and fluid particles respectively.
Density of the wall particle is then set using this pressure
For a two fluid system (boundary, fluid), this equation must be instantiated with boundary as the destination and fluid as the source.
The boundary particle array must additionally define a property \(wij\) for the denominator in Eq. (27) from [Adami2012]. This array sums the kernel terms from the ghost particle to the fluid particle.
Bases: pysph.sph.equation.Equation
Generalized Weakly Compressible Equation of State
This is the generalized Tait’s equation of state and the suggested values in [Adami2013] are \(\mathcal{X} = 0\), \(\gamma=1\) and \(b = 1\).
The reference pressure \(p_0\) is calculated from the artificial sound speed and reference density:
Bases: pysph.sph.equation.Equation
Summation density with volume summation
In addition to the standard summation density, the number density for the particle is also computed. The number density is important for multi-phase flows to define a local particle volume independent of the material density.
For this equation, the destination particle array must define the variable V for particle volume.
Bases: pysph.sph.equation.Equation
Set the inverse volume using mass density
Bases: pysph.sph.equation.Equation
Number density for volume computation
See SummationDensity
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
[Gray2001] | J. P. Gray et al., “SPH elastic dynamics”, Computer Methods in Applied Mechanics and Engineering, 190 (2001), pp 6641 - 6662. |
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Rate of change of stress (2D)
where
Bases: pysph.sph.equation.Equation
Momentum Equation with Artificial Stress
where
Bases: pysph.sph.equation.Equation
Artificial stress to remove tensile instability
The dispersion relations in [Gray2001] are used to determine the different components of \(R\).
Angle of rotation for particle \(a\)
In rotated frame, the new components of the stress tensor are
Components of \(R\) in rotated frame:
Components of \(R\) in original frame:
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Basic equations for Gas-dynamics
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Summation density with iterative solution of the smoothing lengths.
Parameters:
Bases: pysph.sph.equation.Equation
Implementation of the equations used for surface tension modelling, for example in KHI simulations. The references are as under:
- M. Shadloo, M. Yildiz, “Numerical modelling of Kelvin-Helmholtz isntability using smoothed particle hydrodynamics”, IJNME, 2011, 87, pp 988–1006 [SY11]
- Joseph P. Morris “Simulating surface tension with smoothed particle hydrodynamics”, JCP, 2000, 33, pp 333–353 [JM00]
- Adami et al. “A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation”, JCP 2010, 229, pp 5011–5021 [A10]
Bases: pysph.sph.equation.Equation
Gradient of color Eq. (14) in [A10]
where, the average \(\tilde{c}_{ab}\) is defined as
Bases: pysph.sph.equation.Equation
Reproducing divergence approximation Eq. (20) in [A10] to compute the curvature
Bases: pysph.sph.equation.Equation
Acceleration due to surface tension force Eq. (25) in [JM00]:
Note that as per Eq. (17) in [JM00], the un-normalized normal is basically the gradient of the color function. The acceleration term therefore depends on the gradient of the color field.
Bases: pysph.sph.equation.Equation
Gradient of the color function using Eq. (13) of [SY11]:
Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.
Singularities are avoided as per the recommendation by [JM00] (see eqs 20 & 21) using the parameter \(\epsilon\)
Bases: pysph.sph.equation.Equation
Interface curvature using number density. Eq. (15) in [SY11]:
Bases: pysph.sph.equation.Equation
Gradient of the color function using Eq. (17) of [JM00]:
where a smoothed representation of the color is used in the equation. Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.
Singularities are avoided as per the recommendation by [JM00] (see eqs 20 & 21) using the parameter \(\epsilon\)
Bases: pysph.sph.equation.Equation
Gradient of the color function using Eq. (13) of [SY11]:
Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.
Bases: pysph.sph.equation.Equation
Discretized dirac-delta for the SY11 formulation Eq. (14) in [SY11]
This is essentially the same as computing the color gradient, the only difference being that this might be called with a reduced smoothing length.
Note that the normals should be computed using the SY11ColorGradient equation. This function will effectively overwrite the color gradient.
Bases: pysph.sph.equation.Equation
Acceleration due to surface tension force Eq. (7,9) in [SY11]:
where, \(\delta^s\) is the discretized dirac delta function, \(\boldsymbol{n}\) is the interface normal, \(\kappa\) is the discretized interface curvature and \(\sigma\) is the surface tension force constant.
Bases: pysph.sph.equation.Equation
Smoothed color function. Eq. (17) in [JM00]
where, \(c_b^i\) is the color index associated with a particle.
The basic equations for the IISPH formulation of
M. Ihmsen, J. Cornelis, B. Solenthaler, C. Horvath, M. Teschner, “Implicit Incompressible SPH,” IEEE Transactions on Visualization and Computer Graphics, vol. 20, no. 3, pp. 426-435, March 2014. http://dx.doi.org/10.1109/TVCG.2013.105
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
This is important and not really discussed in the original IISPH paper.
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.integrator_step.IntegratorStep
A straightforward and simple integrator to be used for IISPH.
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
The acceleration on the fluid due to a boundary.
Rigid body related equations.
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.integrator_step.IntegratorStep
Fast but inaccurate integrator. Use this for testing
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
The pressure acceleration on the fluid/solid due to a boundary. Implemented from Akinci et al. http://dx.doi.org/10.1145/2185520.2185558
Use this with the fluid as a destination and body as source.
Bases: pysph.sph.integrator_step.IntegratorStep
Bases: pysph.sph.equation.Equation
This is inspired from http://http.developer.nvidia.com/GPUGems3/gpugems3_ch29.html and
BK Mishra’s article on DEM http://dx.doi.org/10.1016/S0301-7516(03)00032-2
A review of computer simulation of tumbling mills by the discrete element method: Part I - contact mechanics
Note that d is a factor multiplied with the “h” of the particle.
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
The viscous acceleration on the fluid/solid due to a boundary. Implemented from Akinci et al. http://dx.doi.org/10.1145/2185520.2185558
Use this with the fluid as a destination and body as source.
Use sympy to perform most of the math and use the resulting formulae to calculate:
inv(I) ( au - w x (I w))
Bases: pysph.sph.equation.Equation
Bases: pysph.sph.equation.Equation
Functions to reduce array data in serial or parallel.
Simply returns the array for the serial case.
Reduce an array given an array and a suitable reduction operation.
Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.
Parameters
- array: numpy.ndarray: Any numpy array (1D).
- op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)
Reduce an array given an array and a suitable reduction operation.
Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.
Parameters
- array: numpy.ndarray: Any numpy array (1D).
- op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)
Reduce an array given an array and a suitable reduction operation.
Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.
Parameters
- array: numpy.ndarray: Any numpy array (1D).
- op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)
Bases: object
A group of equations.
This class provides some support for the code generation for the collection of equations.
Constructor.
When running simulations in parallel, one should typically run the summation density over all particles (both local and remote) in each processor. This is because we must update the pressure/density of the remote neighbors in the current processor. Otherwise the results can be incorrect with the remote particles having an older density. This is also the case for the TaitEOS. In these cases the group that computes the equation should set real to False.