Calculation of Velocity in Post-Processing

Asalih3d · April 2019

Does anyone by chance have an example considering how to calculate this properly?

I am trying to do this on the example with Periodicity in Main at location x = 0.2, y = 0 and z = 0.002, using this code:

set dirout2=%dirout%\measuretool
%measuretool% -dirin %diroutdata% -points points_cyl.txt -savecsv %dirout2%/VelPointBottom2 -vars:-all,+vel -onlytype:-all,+fluid
if not "%ERRORLEVEL%" == "0" goto fail

The result from the measuretool seems a bit weird when compared with my own code:

My own code is basically done in Julia as:

h  = convert(Float32,0.002449) #m, smoothing length for Periodicty, example 2 in Main
x  = convert(Float32,0.2) #y = 0 z = 0 in theory now
#\alphaD from Runout file
αD = 92833.429688 #7/(4*pi*(h/2)^2) #Wendland Quintic 2D
#Abs since a particle on both sides of x should be counted in bool later
q(ra::Float32,rb::Array{Float32,1},h::Float32) = abs.(ra.-rb)/h
#Wendland kernel as in Periodicity
W(q) =  αD*(1 .- q/2).^4 .* (2 .* q .+ 1)
#Reading points data for particles and velocity
pos = readVtkArray("PartFluid_",PostSPH.Points)
vel = readVtkArray("PartFluid_",PostSPH.Vel)
#Initializing x-velocity array, Vab
Vab = []
for i = 1:length(pos)
    qab = q(x,pos[i][:,1],h)
    #Bool is to ensure on values between 0 and 2 are taken into account
    bool = convert(Array{Int},0 .< qab .< 2)
    qabF = q(x,pos[i][:,1],h).*bool
    Wab = W(qabF)
    push!(Vab,sum(Wab.*vel[i][:,1])/sum(Wab))
end

Maybe some one can help me out / show me an example of why I am wrong?

Kind regards

Asalih3d · April 2019

Has no one else tried to calculate manually and verifying everthing is functioning correctly?

Kind regards

Alex · May 2019

The best way to check it is comparing with reference (analytical, experimental) velocity values of an experiment... as we have done in some papers already

And compute velocity far from boundary layer since maybe you are not considering GAP fluid-bound or other boundary effects

Regards

Asalih3d · May 2019

Thanks for your comment. Will try to make the code work for some random point in the domain and see if I can get results, corresponding to DualSPHysics own tool.

Kind regards

Alex · May 2019

Note that we also:
1) normalised the kernel
2) we only compute variables if the sum_wab is higher than a minimum value (0.5 for example)... this helps to avoid computing variables when we have only 2-3 particles in the kernel domain

Alex · May 2019

that is why we use "kclimit" to define the minimum sum_wab (being maxim 1) to normalise, but also to define a dummy value when sum_wab is not big enough (like less than 0.5 for example)
dummy values are "0" for velocity, but can be -999999999 for pressure, etc

Alex

Asalih3d · May 2019

I will try to see if I can find the specific code you used, it must be in one of the files, but thanks so much for your explanation. I am trying to increase the functionality of my post-processing tool, https://github.com/AhmedSalih3d/PostSPH.jl, so it is nice getting some help.

Will let you know if I couldn't find it, in the weekend.

Kind regards

Alex · May 2019

https://github.com/DualSPHysics/DualSPHysics/wiki/11.-Post‐processing#113-analysis-of-numerical-measurements-measuretool:

"Kernel correction is also applied when the summation of the kernel values around the position is higher than a value (-kclimit:) defining a dummy value if the correction is not applied (-kcdummy:). "

Calculation of Velocity in Post-Processing

Comments