Increasing run time as flow behavior index (n) decreases in Power Law fluids

Hi there

I am simulating the Poiseulle flow of Power low fluids by DualSPHysics.

I am using HBP model to stimulate non-Newtonian behavior of Power low fluids by setting the appropriate value of m and n parameters in HBP model. When n ≠ 1 and m = 0 the HBP model reduces to a Power Law model.

I am observing an issue that with the same number of particles, with the same IDP by decreasing flow behavior index (n) the run time increases.

I will be appreciated if someone can explain to me what is the reason behind increasing run time with the same IDP.

Happy new year in advance

Best regards


  • Let us ask @gfourtakas


  • Hi,

    This is to be expected with the fully explicit time integration nature of the code and the viscous CFL condition.

    Let me explain: when the flow behaviour index reduces and the your obtain shear thinning characteristics at the higher end of your strain rate the slope is getting flat i.e. low viscosity. At the other end, low strain rates, your slope is quite sharp (the more you decrease the flow index, the higher the slope), which results in a large viscosity. Plot apparent viscosity vs stress to understand this better.

    As we require to use a very severe CFL condition h^2/ν the timestep suffers. Please use the relaxation parameter (with extreme caution) to try and speed it up a bit :)



  • Dear Dr Fourtakas

    Thank you so much for your response.

    I will try to understand that.

    Merry Christmas

    Best regards

  • Hi

    I also run the Poiseulle flow of Power law fluids on DualSPHysics. I find when I set up n=0.3, it takes too long time. Which parameters I need change? I try to use smaller delt, but failed.

    One quick question, for HBP model, I see documents said when m=10 it means power law fluid, m=100 it means Bingham fluid. It seems diiferernt from your set up. When n ≠ 1 and m = 0 the HBP model reduces to a Power Law model.

    Wish you all the best!


  • Dear Xue

    I suggest you set m=10 if your fluid has a yield stress. Otherwise, setting yield stree=0 or m=0 have the same consequence but they may have different calculation effects.


  • Dear Salehmeiabadi

    Thanks for your replying! I will try your suggestion to see the results.

    Have you ever run cases n<1, if possible can I see the set up of your case?

    Wish you have a good day!

    Best regards,


  • <?xml version="1.0" encoding="UTF-8" ?>




          <gravity x="0.01" y="0" z="0" 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="1" comment="Coefficient to multiply speedsystem" />

          <speedsound value="10" auto="false" comment="Speed of sound to use in the simulation (by default speedofsound=coefsound*speedsystem)" />

          <coefh value="1.0" comment="Coefficient to calculate the smoothing length (h=coefh*sqrt(3*dp^2) in 3D)" />

          <cflnumber value="0.1" comment="Coefficient to multiply dt" />


        <mkconfig boundcount="240" fluidcount="9"/>


          <definition dp="0.00002">

            <pointmin x="-0.00005" y="0" z="-0.0015" />

            <pointmax x="0.01" y="0" z="0.0015" />




              <setshapemode>dp | bound</setshapemode>

              <setdrawmode mode="full" />

    <setmkfluid mk="1" />



                <point x="0" y="-0.05" z="-0.0005" />

                <size x="0.0015" y="0.1" z="0.001" />


              <setmkbound mk="1" />



                <point x="0" y="-0.05" z="-0.0005" />

                <size x="0.0015" y="0.1" z="0.001" />


              <shapeout file="Box" />







          <nnphases> %Defines non-newtonian phases parameters

            <phase mkfluid="1">

              <rhop value="1000" comment="Density of the phase" />           

              <visco value="1e-6" comment="Kinematic viscosity (or consistency index) value for phase (m2/s)" />

              <tau_yield value="0" comment="User defined maximum 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="0.95" 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" />





          <parameter key="SavePosDouble" value="0" comment="Saves particle position using double precision (default=0)" />

          <parameter key="StepAlgorithm" value="1" 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="2" 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)" />

          <parameter key="TimeMax" value="0.8" comment="Time of simulation" units_comment="seconds" />

          <parameter key="TimeOut" value="0.05" comment="Time out data" units_comment="seconds" />

          <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" />

    <parameter key="XPeriodicIncZ" value="0" comment="Increase of Z with periodic BC" />

          <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%" />





  • I had the same problem. I still can't understand why n < 1, it is so slow. Who can tell me why?

  • Since when shear thinning occur, viscosity will be smaller. So the time step should be larger?

  • @gfourtakas can help here

  • Hi TianwenDong and all,

    Let me clear this up for you.

    For a shear thinning material at low shear (i.e. before it becomes thinning) there is quite a large gradient that is in other words viscosity tends to high values

    I have been overcautious on this matter with this implementation and since we are using a fully explicit formulation the viscous CFL condition is severe.

    I do not recommend using such low values of n i.e., n<0.3 and I have never tested below n = 0.5-0.6 for the problems I was looking at.

    Thank you


Sign In or Register to comment.