Kelvin-Helmholtz instabilities

idaapm

Hi,

I'm trying to reproduce the Kelvin-Helmholtz instabilities based on Thorpe's experiment (1969) with DualSPHysics.

This is my XML file:
- - -
case app="GenCase v2.03 (27-09-2013)" date="01-02-2016">
casedef>
constantsdef>
lattice bound="2" fluid="1" />
gravity x="0" y="0" z="-9.81" />
cflnumber value="0.1" />
hswl value="0" auto="true" />
coefsound value="20" />
coefficient value="0.92" />
gamma value="7" />
rhop0 value="1000" />
eps value="0.0" />
constantsdef>
mkconfig boundcount="240" fluidcount="10">
geometry>
definition dp="0.001">
pointmin x="-10" y="0.015" z="-10" />
pointmax x="10" y="0.015" z="10" />
definition>
commands>
mainlist>
setshapemode>dp | bound |fluid
setdrawmode mode="full" />
setmkfluid mk="0"/>
matrixreset/>
rotate x="0" y="1" z="0" ang="4.13"/>
drawbox>
boxfill>solid
point x="0.000" y="0.000" z="0.000" />
size x="1.830" y="0.030" z="0.015" />
drawbox>
setmkfluid mk="1"/>
matrixreset/>
rotate x="0" y="1" z="0" ang="4.13"/>
drawbox>
boxfill>solid
point x="0.000" y="0.000" z="0.015" />
size x="1.830" y="0.030" z="0.015" />
drawbox>
setdrawmode mode="full" />
setmkbound mk="1" />
matrixreset/>
rotate x="0" y="1" z="0" ang="4.13"/>
drawbox>
boxfill>all
point x="0.000" y="0.000" z="0.000" />
size x="1.830" y="0.0300" z="0.030" />
drawbox>
shapeout file=""/>
mainlist>
commands>
geometry>
casedef>
execution>
parameters>
parameter key="StepAlgorithm" value="2" comment="Step Algorithm 1:Verlet, 2:Symplectic (def=1)" />
parameter key="VerletSteps" value="40" comment="Verlet only: Number of steps to apply Eulerian equations (def=40)" />
parameter key="Kernel" value="2" comment="Interaction Kernel 1:Cubic Spline, 2:Wendland (def=1)" />
parameter key="ViscoTreatment" value="2" comment="Viscosity Formulation 1:Artificial, 2:Laminar+SPS (def=1)" />
parameter key="Visco" value="0.0010" comment="Viscosity value" />
parameter key="ShepardSteps" value="0" comment="Number of steps to apply Shepard density filter, 0=non applied (def=0)" />
parameter key="DeltaSPH" value="0.1" comment="DeltaSPH value, 0.1 is the typical value, with 0 disabled (def=0)" />
parameter key="DtIni" value="0.000000001" comment="Initial time step" />
parameter key="DtMin" value="0.0000000001" comment="Minimum time step (def=0.00001)" />
parameter key="TimeMax" value="20.0" comment="Time of simulation" />
parameter key="TimeOut" value="0.01" comment="Time out data" />
parameter key="IncZ" value="5" comment="Increase of Z+" />
parameter key="PartsOutMax" value="1" comment="Proportion of fluid particles out allowed (def=1)" />
parameters>
execution>
case>
- - -

I've also changed both JSphCpu and JSphCpuSingle files in order to use different viscosity and density values for each fluid:
- - -
JSphCpu.cpp @ Laminar+SPS viscosity
if(Idp[p]>=Npb && Idp[p]=Npb+25294)Visco=0.0015f;

JSphCpuSingle.cpp @ Configuration of the current domain
for(unsigned p=0;p=Npb && p=Npb+25294)Rho[p]=780.f;
}
- - -

The results are not satisfactory. I should expect both fluids to be moving in the different directions (at least in the beginning of the simulation) and the interface between them to stabilize horizontally. I'm not getting neither.

I've also noticed that the liquid seems to be leaking and that there's an "empty" space that appears in one of the corners. Where does the "empty" space come from? And why is my fluid leaking?

Is there a way of solving this problem?

Thanks.

Kind regards, ida a.


This is my Run.out file:
- - -
DualSPHysics v3 (17-12-2013)
=============================
Threads by host for parallel execution: 4
[Initialising JSphCpuSingle v3.00 03-02-2016 22:01:35]
**Case configuration is loaded
Loading initial state of particles...
Loaded particles: 61714 (61714+0)
MapPos(border)=(-0.002251,0.014999,-0.000251)-(1.825251,0.015001,0.162251)
MapPos=(-0.002251,0.014999,-0.000251)-(1.825251,0.015001,0.974764)
**Initial state of particles is loaded
**2D-Simulation parameters:
CaseName="CasePeriodicity"
RunName="CasePeriodicity"
SvTimers=True
StepAlgorithm="Symplectic"
Kernel="Wendland"
Viscosity="LaminarSPS"
Visco=0.001000
ShepardSteps=0
DeltaSph="DBCExt"
DeltaSphValue=0.100000
CaseNp=61714
CaseNbound=8932
CaseNfixed=8932
CaseNmoving=0
CaseNfloat=0
PeriodicActive=0
Dx=0.001000
H=0.001301
CteB=16256.571289
Gamma=7.000000
Rhop0=1000.000000
Eps=0.000000
Cs0=10.667521
CFLnumber=0.100000
DtIni=0.000000
DtMin=0.000000
MassFluid=0.001000
MassBound=0.000500
WendlandCte.awen=329040.593750
WendlandCte.bwen=-1264584576.000000
SpsSmag=0.000000
SpsBlin=0.000000
TimeMax=20.000000
TimePart=0.010000
Gravity=(0.000000,0.000000,-9.810000)
PartOutMax=52782
RhopOut=True
RhopOutMin=700.000000
RhopOutMax=1300.000000
CellOrder="XYZ"
CellMode="2H"
Hdiv=1
MapCells=(703,1,375)
RunMode="OpenMP(Threads:4, mode:Static)"
Allocated memory in CPU: 24283896 (23.16 MB)
Part_0000 61714 particles successfully stored


Particles out: 4 (total: 4)

Particles out: 4 (total: 8)

Particles out: 1 (total: 9)

Particles out: 5 (total: 14)

Particles out: 1 (total: 15)

Particles out: 4 (total: 19)

Particles out: 1 (total: 20)

Particles out: 1 (total: 21)

[Simulation finished 07-02-2016 04:44:55]
DTs adjusted to DtMin.............: 0
Excluded particles................: 21
Total Runtime.....................: 283400.031250 sec.
Simulation Runtime................: 283400.031250 sec.
Time per second of simulation.....: 14169.997070 sec.
Steps per second..................: 7.311636
Steps of simulation...............: 2072118
PART files........................: 2001
Maximum number of particles.......: 61714
Maximum number of cells...........: 44289
CPU Memory........................: 26509872 (25.28 MB)

Alex

If you visit:

You can read:

The RT instability occurs when a heavy fluid (1800 kg/m^3) lies on top of a light fluid (1000 kg/m^3) in the attendance of a gravitational field. This instability phenomenon has a relevant importance in many industrial processes. Intial perturbation, z=0.15sin(2ᴨx).

DualSPHysics includes here:
- An inter-particle interaction force included to improve the interface treatment. This approach was proposed by Tartakovsky & Meakin (2005) to simulate a surface tension effect and also wettability effects. The inter-particle interaction force is added in the equation of motion.
- Variable smoothing length was proposed for Sigalotti & Lopez (2008) as a mean to stabilize and properly describe the interface of fluids.
- Shifting algorithm proposed by Xu et al. (2009) and it is used to keep a better distribution of particles. The magnitude and direction of the position shift is governed by Fick’s law, which slightly moves particles from high concentration regions to regions with less number of particles.

Implementation by Carlos Enrique Alvarado Rodríguez, Universidad de Guanajuato (Mexico).

So I would suggest you to email Carlos in iqcarlosug@gmail.com since he is preparing a paper where details will be presented.

Regards

idaapm

Many thanks Alex! I'll contact him for further information about my issue.

idaapm

Hi again,

I've been implementing changes in the source code (version 4) to be able to simulate a fluid-fluid flow.

I'm testing a 2D square cavity (1x1 m²) partially filled with water (at the bottom, h=0.4 m) and oil (at the top, h=0.4 m). I've set dp= 0.01 m. I've modified the gravity vector so that I have a rotated cavity: <gravity x="0.71" y="0" z="-9.78". My fluid properties are the following: water density=1000 kg/m³ viscosity=0.001 kg/m.s and oil density=780 kg/m³ viscosity=0.0015 kg/m.s .

Depending on the particle id I set different properties for the fluid phases. The changes account for:

@JSph.cpp: RhopZero + MassFluid + MassBound

@JSphCpu.cpp: visco + RhopZero (for the EoS)

@JSphCpuSingle.cpp: Velrhopc[p].w (at ConfigDomain)

I've tested different options, but I don't get any imrpovments. I get an initial explosion in both phases and some particles are excluded from my initial domain. The bottom fluid interface doesn't behave as expected (it should be parallel to the top one).

Could you help me understand what is happening with my simulations? Do you have any suggestion or tips to solve this problem? Can the explosion be related with the distance between boundary and fluid particles? Could it also be related with the implementation of both MassFluid and MassBound?

Kind regards, ida a.

Alex

As mentioned before you have to contact Carlos Enrique Alvarado Rodríguez, Universidad de Guanajuato (Mexico) in iqcarlosug@gmail.com.

He is the main author of the work:



A journal paper is being produced but you can also ask him for the proceedings:

Alvarado-Rodriguez, C.E., Barreiro, A., Dominguez, J.M., Crespo, M., Gómez-Gesteira, M., Klapp, J., 2015. Simulation of dispersion in a porous media with multiphase Smoothed Particle Hydrodinamics. In: Proceedings of the 10th SPHERIC International Workshop, Parma, Italy.

Regards

idaapm

Many thanks, Alex. Will get in contact with Carlos again.

Alex

You can also email us (dualsphysics@gmail.com) and we will send to you that proceedings paper

Regards

idaapm

Will do. Thanks!

Ranit

Hello,

I am trying to simulate KH Instability, but I am getting an error that Execution with Periodic Boundary Conditions is not available. I have based my set-up on CaseTwoPhases.

Thanks in Advance

Alex

That is not an error... It is information for the user!!
The version 3.2 of the executable that includes the multiphase implementation DO NOT include periodic boundaries.
Source code will be released soon. If periodicity is not implemented, it can be included by the user in the code.

Regards

Ranit

Thanks Alex.

Also, I am trying to have an initial disturbance (of displacement) at the interface by importing initial geometry of fluids from an external .stl file.

Is there a better way to create initial disturbance at the interface within DualSPHysics?

Regards
Ranit