Help with Sand Simulation Behavior in DualSPHysics

[Murat Gokoglu]

Howdy All,

I am working on a simulation involving both water and sand. I've set up my model, but I am having trouble with the sand. When I run the simulation, the sand doesn't behave like sand — it seems to act more like a fluid. Below is an excerpt of the code I am using for defining the sand phase:

<?xml version="1.0" encoding="UTF-8" ?>
<case app="GenCase v5.0.164 (21-03-2020)" date="26-03-2020 02:22:19">
    <casedef>
        <constantsdef>
            <gravity x="0" y="0" z="-9.81" comment="Gravitational acceleration" units_comment="m/s^2" />
            <rhop0 value="1000" comment="Reference density of the fluid" units_comment="kg/m^3" />
            <hswl value="0" auto="true" comment="Maximum still water level to calculate speedofsound using coefsound" units_comment="metres (m)" />
            <gamma value="7" comment="Polytropic constant for water used in the state equation" />
            <speedsystem value="0" auto="true" comment="Maximum system speed (by default the dam-break propagation is used)" />
            <coefsound value="20" comment="Coefficient to multiply speedsystem" />
            <speedsound value="0" auto="true" comment="Speed of sound to use in the simulation (by default speedofsound=coefsound*speedsystem)" />
            <coefh value="1.5" comment="Coefficient to calculate the smoothing length (h=coefh*sqrt(3*dp^2) in 3D)" />
            <cflnumber value="0.2" comment="Coefficient to multiply dt" />
        </constantsdef>
        <mkconfig boundcount="230" fluidcount="9">
            <!--mkorientfluid mk="0" orient="Xyz" />
            <mkorientfluid mk="1" orient="Xyz" /-->
        </mkconfig>
        <geometry>
            <!--DEFINITION OF DOMAIN WHERE PARTICLES WILL BE CREATED -->
            <definition dp="0.05" comment="Initial inter-particle distance" units_comment="metres (m)">
                <pointref x="0.0" y="0.0" z="0.0" />
                <pointmin x="0.0" y="0.0" z="0.0" />
                <pointmax x="25.0" y="0.0" z="10.0" />
            </definition>
            <commands>
                <mainlist>
                <setshapemode>actual | dp | bound</setshapemode>
                <setmkbound mk="2"/>
                <setdrawmode mode="full"/>
                <drawbox objname="Piston">
                    <boxfill>solid</boxfill>
                    <point x="0.1" y="-2.5" z="0.0" />
                    <size x="0.01" y="5.0" z="3.0" />
                </drawbox>
                    <!--CREATION OF BOUNDARY PARTICLES (WALLS OF TANK) -->
                    <setmkbound mk="0"/>
<setdrawmode mode="full"/>
<drawbox objname="Bottom">
<boxfill>solid</boxfill>
<point x="0.0" y="-2.5" z="0.0" />
<size x="20.0" y="5.0" z="0.01" />
</drawbox>
<move x="0.0" y="-2.5" z="0.0" />
<setmkfluid mk="8"/>
<fillbox x="5.0" y="2.5" z="1.0" objname="FillBox">
                    <modefill>void</modefill>
<point x="0" y="0" z="0" />
<size x="8.0" y="5.0" z="1.5" />
</fillbox>
<matrixreset />
<shapeout file="" />
<!-- making a polygon for sand -->
<!--setshapemode>dp | bound</setshapemode-->
<setmkfluid mk="0" />
<!--fillprism x="405" y="0" z="13"-->
<!--setdrawmode mode="solid" /-->
<!--drawextrude closed="true"-->
<setmkfluid mk="0" /> <!-- Marker for Sand -->
<setmkfluid mk="0" /> <!-- Marker for Sand -->
<drawprism>
<!-- Bottom surface of the sand dune -->
<point x="10.0" y="0" z="0.09" />
<point x="12.0" y="0" z="0.09" />
<point x="14.0" y="0" z="0.09" />

<!-- Top surface of the sand dune -->
<point x="10.0" y="0" z="0.5" />  <!-- Starting slope of the dune -->
<point x="12.0" y="0" z="2" />  <!-- Peak of the dune -->
<point x="14.0" y="0" z="0.5" />  <!-- End of the dune -->
</drawprism>


                    <!--CREATION OF BOUNDARY PARTICLES (STRUCTURE FROM STL => stl IS 8 M BELOW MWL) -->
                   
                    <!--CREATION OF FLUID PARTICLES (FILLBOX WITH WATER=> stl IS PLACED 8 M BELOW MWL) -->
  <!--setmkfluid mk="8" />
                    <fillbox x="1.5" y="-0.2" z="8.0">
                        <modefill>void</modefill>
                        <point x="1.8" y="0" z="1" />
                        <size x="15" y="0.37" z="9.0" />
                    </fillbox-->
                </mainlist>
            </commands>
        </geometry>
        <motion>
            <objreal ref="2">
                <begin mov="1" start="0"/>
                <mvnull id="1" />
            </objreal>
        </motion>
    </casedef>
    <execution>
<special>
        <initialize>
        </initialize>
            <wavepaddles>
                <piston>
                    <mkbound value="2" comment="Mk-Bound of selected particles" />
                    <start value="0" comment="Start time (default=0)" />
                    <duration value="0.8" comment="Movement duration, Zero is the end of simulation (default=0)" />
                    <depth value="0.8" comment="Water depth (default=0)" />
                    <pistondir x="1.0" y="0.0" z="0.0" comment="Movement direction (default=(1,0,0))" />
                    <waveorder value="2" comment="Order wave generation 1:1st order, 2:2nd order (default=1)" />
                    <waveheight value="0.3" comment="Wave height" />
                    <waveperiod value="1.0" comment="Wave period" />
                    <gainstroke value="1.0" comment="Gain factor to amplify/reduce the paddle stroke (default=1)" />
                    <phase value="0.0" comment="Initial wave phase in function of PI (default=0)" />
                    <ramp value="0.0" comment="Periods of ramp (default=0)" />
                    <savemotion periods="24" periodsteps="20" xpos="2.0" zpos="-0.15" comment="Saves motion data. xpos and zpos are optional. zpos=-depth of the measuring point" />
                </piston>
            </wavepaddles>
            <nnphases> %Defines non-newtonian phases parameters
                <phase mkfluid="8">
                   <rhop value="1000" comment="Density of the phase" />                    
                    <visco value="0.001" comment="Kinematic viscosity (or consistency index) value for phase (m2/s)" />
<tau_yield value="0.0005" comment="Specific yield stress of phase (Pa m3/kg)" />
                    <HBP_m value="0" comment="Use 0 to reduce Newtonian liquid, order of 10 for power law and order of 100 for Bingham (sec) " />
                    <HBP_n value="1" comment="Use 1 to reduce to Newtonian, HBP_n<1 for shear thinning HBP_n>1 for shear thickenning " />    
<phasetype value="0" comment="Non-Newtonian=0 only option in v5.0" />                
                    <!--rhop value="1500" comment="Density of the phase" />
                    <visco value="0.01" comment="Kinematic viscosity (or consistency index) value for phase (m2/s)" />
                    <tau_yield value="0.0001" comment="Specific yield stress of phase (Pa m3/kg)" />
<tau_max value="0.00001" comment="User defined maximum specific yield stress of phase (Pa m3/kg)" />
<Bi_multi value="10.0" comment="tau_max multiplier for use with Bingham model or bi-viscosity model(tau_bi=tau_max*Bi_multi)" />
                    <HBP_m value="10" comment="Use 0 to reduce Newtonian liquid, order of 10 for power law and order of 100 for Bingham (sec) " />
                    <HBP_n value="1.5" comment="Use 1 to reduce to Newtonian, <1 for shear thinning >1 for shear thickenning " />                    
<phasetype value="0" comment="Non-Newtonian=0 only option in v5.0" /-->
               </phase>  
                <phase mkfluid="0">
                   <rhop value="2690" comment="Density of the phase" />                    
                    <visco value="50" comment="Kinematic viscosity (or consistency index) value for phase (m2/s)" />
                    <tau_yield value="0.0001" comment="Specific yield stress of phase (Pa m3/kg)" />
<tau_max value="10" comment="User defined maximum specific yield stress of phase (Pa m3/kg)" />
<Bi_multi value="10.0" comment="tau_max multiplier for use with Bingham model or bi-viscosity model(tau_bi=tau_max*Bi_multi)" />
                    <HBP_m value="100" comment="Use 0 to reduce Newtonian liquid, order of 10 for power law and order of 100 for Bingham (sec) " />
                    <HBP_n value="1.5" comment="Use 1 to reduce to Newtonian, <1 for shear thinning >1 for shear thickenning " />                    
<phasetype value="1" comment="Non-Newtonian=0 only option in v5.0" />
                </phase>  
            </nnphases>
        </special>
        <parameters>
            <parameter key="SavePosDouble" value="0" comment="Saves particle position using double precision (default=0)" />
            <parameter key="StepAlgorithm" value="2" comment="Step Algorithm 1:Verlet, 2:Symplectic (default=1)" />
            <parameter key="VerletSteps" value="40" comment="Verlet only: Number of steps to apply Euler timestepping (default=40)" />
            <parameter key="Kernel" value="2" comment="Interaction Kernel 1:Cubic Spline, 2:Wendland (default=2)" />
%Choice of reology treatment, velocity gradient calculation and viscosity treatment  
<parameter key="RheologyTreatment" value="2" comment="Reology formulation 1:Single-phase classic, 2: Single and multi-phase" />
            <parameter key="VelocityGradientType" value="1" comment="Velocity gradient formulation 1:FDA, 2:SPH" />
            <parameter key="ViscoTreatment" value="1" comment="Viscosity formulation 1:Artificial, 2:Laminar+SPS, 3:Constitutive  eq." />
            %Wall boundary viscosity or/and artificial viscosity if ViscoTreatment is 1:Artificial
<parameter key="Visco" value="0.05" comment="Viscosity value" /> % Note alpha can depend on the resolution when using artificial viscosity
            <parameter key="ViscoBoundFactor" value="1" comment="Multiply viscosity value with boundary (default=1)" />
            <parameter key="DensityDT" value="3" comment="Density Diffusion Term 0:None, 1:Molteni, 2:Fourtakas, 3:Fourtakas(full) (default=0)" />
            <parameter key="DensityDTvalue" value="0.1" comment="DDT value (default=0.1)" />
            <parameter key="Shifting" value="3" comment="Shifting mode 0:None, 1:Ignore bound, 2:Ignore fixed, 3:Full (default=0)" />
            <parameter key="ShiftCoef" value="-10" comment="Coefficient for shifting computation (default=-2)" />
<parameter key="ShiftTFS" value="1.5" comment="Threshold to detect free surface. Typically 1.5 for 2D and 2.75 for 3D (default=0)" />
            <parameter key="RigidAlgorithm" value="1" comment="Rigid Algorithm 0:collision-free, 1:SPH, 2:DEM, 3:Chrono (default=1)" />
            <parameter key="FtPause" value="0.0" comment="Time to freeze the floatings at simulation start (warmup) (default=0)" units_comment="seconds" />
            <parameter key="CoefDtMin" value="0.05" comment="Coefficient to calculate minimum time step dtmin=coefdtmin*h/speedsound (default=0.05)" />
            <!--parameter key="RelaxationDt" value="0.2" comment="Relaxation parameter for the viscous time step restricition(default=0.2)" /-->
<parameter key="DtIni" value="0" comment="Initial time step. Use 0 to defult use (default=h/speedsound)" units_comment="seconds" />
            <parameter key="DtMin" value="0" comment="Minimum time step. Use 0 to defult use (default=coefdtmin*h/speedsound)" units_comment="seconds" />
            <parameter key="DtFixed" value="0" comment="Fixed Dt value. Use 0 to disable (default=disabled)" units_comment="seconds" />
            <parameter key="DtFixedFile" value="NONE" comment="Dt values are loaded from file. Use NONE to disable (default=disabled)" units_comment="milliseconds (ms)" />
            <parameter key="DtAllParticles" value="0" comment="Velocity of particles used to calculate DT. 1:All, 0:Only fluid/floating (default=0)" />
             <!-- time   TIME IS IN SECONDS, ParaView has more steps because of TimeOut-->
            <!--parameter key="TimeMax" value="100" comment="Time of simulation" /-->
            <parameter key="TimeMax" value="10" comment="Time of simulation" />
            <parameter key="TimeOut" value="0.1" comment="Time between output files" />
            <parameter key="PartsOutMax" value="1" comment="%/100 of fluid particles allowed to be excluded from domain (default=1)" units_comment="decimal" />
            <parameter key="RhopOutMin" value="500" comment="Minimum rhop valid (default=700)" units_comment="kg/m^3" />
            <parameter key="RhopOutMax" value="3000" comment="Maximum rhop valid (default=1300)" units_comment="kg/m^3" />
            <simulationdomain comment="Defines domain of simulation (default=Uses minimun and maximum position of the generated particles)">
                <posmin x="default" y="default" z="default" comment="e.g.: x=0.5, y=default-1, z=default-10%" />
                <posmax x="default" y="default" z="default + 50%" />
            </simulationdomain>
</parameters>
    </execution>
</case>

Thank you

[Bonaventura TAGLIAFIERRO]

Dear DualSPHysics user,

Thank you for posting regarding your concerns on fluid-to-soil interaction problem. Note, however, that this group is mainly focused on supporting Project Chrono users.

I suggest you report your DualSPHysics related issue here: https://forums.dual.sphysics.org/discussions

If you need any further assistance, please don't hesitate to ask.

Best regards,

Bonaventura