DDT Fourtakas - hydrostatic component

edited December 2021 in DualSPHysics v5.0

Hi,

I'm a bit confused with Fourtakas DT. I know I need to obtain hydrostatic part of density from hydrostatic pressure difference p_ij = rho*g*z_ij, I have to use EOS

for that and then substract rho0. But presented equations contain some weird fractions and square roots.

DualSPHysics wiki says:


struggled with finding of the original reference paper [1], but I looked into [2] instead and there is:

I still haven't understood this formula, so I "asked" DualSPHysics source code:

EDIT: I should add: DDTgz=float(double(RhopZero)*double(fabs(Gravity.z))/double(CteB));

this seems ok to me. But I would like to ask how the previous two equations should work, what I am doing wrong or if there is any notation misinterpretation (or typo?). 

Thank you for your answers.

REFS:

[1]: Fourtakas G, Domínguez JM, Vacondio R, Rogers BD. 2019. Local Uniform Stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models. Computers & Fluids, 190: 346-361. doi.org/10.1016/j.compfluid.2019.06.009.


[2]: Fourtakas G, Domínguez JM, Vacondio R, Rogers BD. 2019. Local Uniform Stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models. Computers & Fluids, 190: 346-361. doi.org/10.1016/j.compfluid.2019.06.009.

Comments

  • If you want to obtain the hydrostatic pressure, I have had success with:

    P_s = rho*g*z_i

    And you can convince your self of this by modelling a still standing fluid with a fixed height.

    Perhaps I did not understand your question fully though.

    Kind regards

  • Perhaps @gfourtakas can help here.

  • Hello Tomas!

    I want first and foremost to apologize for my very unsatisfying answer a year ago. I hope you were able to solve your issue, else I believe I found the proper solution now.

    I agree with you fully that the equation should be written as done in the DualSPHysics code and NOT as in the wiki, which is very subpar in this section. We can easily make sense of this since we know that rho0*g is always roughly equal to 10000. Then if we had to say (rho0*g*drz + 1)/Cb, then the moment drz is "just" slightly negative, i.e. -0.01, we get a negative square root. This has to be an error in the paper.

    Unfortunately there was another mistake in the original paper by Fourtakas, which is in:


    To get the correct "direction" on the density transfer, it should have matched the dynamic density equation. To be honest with you I am not sure how they derived this equation, it seems to me to be a "good idea" taking the principle of artificial viscosity and moving it into density diffusion. The reason I found this inverted ij mistake was due to by chance stumbling on this great paper:

    https://arxiv.org/abs/2110.10076

    Stability and accuracy of the weakly compressible SPH with particle regularization techniques

    Mojtaba JandaghianHerman Musumari SiabenAhmad Shakibaeinia


    Using this formulation I was able to get it to work, see result at the end. The code I used to do it, just as a pointer:

    And the difference in results shown below:


    My own simulation of a dam break. I have particle clumping issues to the left ignore that for now. I am using DDT=2 equivivalent code. As you can see this formulation works and really smooths out the pressure/density field in such a beautiful way!

    @gfourtakas first I want to thank you and the others who have looked into methods of stabilizing the density field in SPH, especially weakly-compressible SPH. This idea of density diffusion has allowed SPH to step up a level in regards to what is possible to simulate and get "clear" results. If I may, I would politely ask you if you would be willing to chime in about what equations are the "true" equations and how to get to them?

    Kind regards

  • @TomasCTU

    Some of my notes to derive the expression and prove it is correct:


    @gfourtakas feel free to comment if I done it wrong or some wrong explanations :)

    Kind regards

Sign In or Register to comment.