I have implemented a constitutive equation of elastic materials (Hooke's law) in my 3D weakly compressible SPH solver based on [1]. The coding seems to be correct. To verify the implementation I started to check the solver through simple test cases. One of them is a rotating planar object in zero gravity (initial conditions contain the velocity field of rigidbody-like rotation).

The problem:

The offdiagonal members of the Cauchy-stress tensor (concerning each particle) should remain zero, but they start to rise as the body deviates from its initial configuration. The larger the deviation the larger the shear-stress in the particles. This results in an oscillatory motion of the body instead of pure rotation. In [1] very large deflections occur, however in my solver these are obstructed by the spurious shears.


What is the reason of this behaviour? Due to the Jaumann-rate the stress-derivative should be objective. Is it true? Are there any non-trivial steps in the calculations of [1] that are not introduced?


[1] J.P. Gray, J.J. Monaghan, R.P. Swift. SPH elastic dynamics, 2001. URL:http://www.sciencedirect.com/science/article/pii/S0045782501002547


You may want to check whether angular momentum is conserved for each particle throughout the simulation.

The artificial viscosity introduced by the model described in that paper, while useful to dampen spurious modes from tensile instability, may have altered the conservation properties of the cluster of SPH particles representing your object.

  • $\begingroup$ Thank you for the contribution. The rotation is obstructed by the Cauchy stress tensor even without artificial viscosity. The offdiagonal elements should remain close to zero in case of a rigid body rotation, but they start to increase as the body rotates forcing the body to stop and rotate back to its original state. $\endgroup$
    – BalazsToth
    Sep 13 '20 at 21:46
  • $\begingroup$ Perhaps try to keep a total Lagrangian kernel approach: determine the neighbours of each particle at the beginning of the simulation, and never update them at any step afterwards. You should also check that for each particle the normalised sum of the kernel shape function over the neighbours is 1, and that the sum of kernel first derivative over neighbours for each particle is instead 0. If they are not conserved, you should apply kernel correction coefficients, I can list papers on the subject. $\endgroup$
    – Joger
    Sep 13 '20 at 21:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.