Stability of pressures in gas


First of all, apologizes from a newbie if I have passed over something obvious.

I would like to share and discuss some behaviour of gas phases I have observed in DualSPHysics4.0_LiquidGas distributed within the DualSPHysics 4.2 package.

It is mostly related with pressure. During simulation of an Oscillating Water Column energy device I suspected that pressure results took time to stabilize, and when it finally does the results seem incorrect.

Absolute pressure in the gas is defined when cs0, gamma and rhop0 is defined as speed of sound can be expressed the following way:

cs0 = (gamma*p0/rhop0)^0.5

So, absolute pressure can be deduced:

p0 = cs0^2*rhop0/gamma

DualSPHysics formulation for pressure is described in chapters 3.4 and 3.13 of "DualSPHysics_v4.2_GUIDE.pdf":

p = (cs0^2*rhop0/gamma) *((rhop/rhop0)^gamma-1)

As cs0^2*rhop0/gamma = p0 this can be written in the following manner:

(p0+p)/p0 = (rhop/rhop0)^gamma

That is the equation for isentropic processes, where p is not the absolute pressure, but the pressure over the main ambient pressure p0, or what is called in engineering as "gauge pressure" pg .

After some intermediate tests I reached to a very simple simulation to isolate ideal gas behaviour in DualSPHysics4.0_LiquidGas that is defined in a later posted XML definition case and can be run by anyone.
The bounday is a column shaped closed chamber 1 m heigh. It is filled with gas particles whose thermodynamic state is defined by the following variables:

rhop0 value="1.204"
gamma value="1.4"
cs0 value="343.28"

Those variables describe air at normal conditions (approx: 20ºC/273K, 100 kPa).

We should expect that after running this case, the air in the upper of column chamber should havean gauge pressure of 0 kPa. The lower gauge pressure should be 11.8 Pa higher due to the weight of the air column (1.2*9.81), and a gauge pressure gradient should be observed in the intermediate layers. And this results should be static, as no motion nor transient behaviour is defined. General pressure p0 is not postprocessed by the code, and so no field of absolute pressure is available, though it can be externally postprocessed.

When this case is run, we can see in a first output (first image) that the code has defined an initial state where upper gauge pressure is 0 and a correct 11.8 Pa gradient is defined.

In a second output (second image), a general increase of gauge pressure is observed: A general gauge pressure of 980 Pa is present with an additional gradient of 11.8 Pa.

The following time steps, main gauge pressure continue a more gradual increase until time step 960 approx (corresponding to 96 seconds, third image), where general gauge pressure stabilizes at 1644 Pa and no apparent gradient. Density results look correct.

1) Why this gradual increase in general gauge pressure?

2) Why this time of 96 seconds for this change in gauge pressure?
Changes of pressure should be governed by speed of sound and would take 1/343 = 0.003 s to be transmitted in a 1 m chamber.

3) Why no pressure gradient remains?

Anyone can help?


P.D. I attach the XML definition case in a following post.


  • Find attached the compressed XML case definition file.
  • I post again the XML file as it is missing:

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

    <case date="06/11/2018 16:42:25">



               <lattice bound="1" fluid="1"/>

               <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" /> %Not used for liquid-gas flows. Set in the properties section for each material.

               <hswl value="0" auto="true" comment="Maximum still water level to calculate speedofsound using coefsound" units_comment="metres (m)" /> %Not used for liquid-gas flows.

               <gamma value="7" comment="Polytropic constant for water used in the state equation" />   %Not used for liquid-gas flows. Set in the properties section for each material.

               <speedsystem value="0" auto="true" comment="Maximum system speed (by default the dam-break propagation is used)" /> %Not used for liquid-gas flows.

               <coefsound value="20" comment="Coefficient to multiply speedsystem" /> %Not used for liquid-gas flows.

               <speedsound value="0" auto="true" comment="Speed of sound to use in the simulation (by default speedofsound=coefsound*speedsystem)" /> %Not used for liquid-gas flows. Set in the properties section for each material.

               <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="20" fluidcount="2"/>


               <definition dp="0.01">

                   <pointmin x="-0.1" y="0" z="-0.1" />

                   <pointmax x="0.2" y="0" z="1.1" />




                       <setshapemode>real | dp | fluid | bound | void</setshapemode>

                       <setdrawmode mode="full" />

                       <setmkfluid mk="1" name="Air" />



                           <point x="0.010" y="-0.01" z="0.010" />

                           <size x="0.080" y="0.01" z="0.980" />


                       <shapeout file="Air" />

                       <setdrawmode mode="face" />

                       <setmkbound mk="0" name="Wall" />


                           <boxfill>left | top | bottom | right </boxfill>

                           <point x="0.000" y="-0.01" z="0.000" />

                           <size x="0.100" y="0.01" z="1.000" />


                       <shapeout file="Wall" />




       <properties> %Material properties for liquid-gas flows are set here.


             <link mkbound="0,1,2" property="water"/>

             <link mkfluid="1" property="air" />


          <property name="water">

             <rhop0 value="1000"/>

             <gamma value="7"/>

             <cs0 value="24.14"/>

             <visco value="0.03"/>

             <phasetype value="0"/>


          <property name="air">

             <rhop0 value="1.204"/>

             <gamma value="1.4"/>

             <cs0 value="343.28"/>

             <visco value="0.03"/>

             <phasetype value="1"/>






               <parameter key="PosDouble" value="1" comment="Precision in particle interaction 0:Simple, 1:Double, 2:Uses and saves double (default=0)" />

          <parameter key="FlowType" value="2" comment="Select Single-phase (1) or liquid-gas flow(2)" />         

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

               <parameter key="ViscoTreatment" value="1" comment="Viscosity formulation 1:Artificial, 2:Laminar+SPS (default=1)" />

               <parameter key="Visco" value="0.0" comment="Viscosity value" /> % Note it is not used for multi-phase flows. The "visco" variable in the properties section will be used instead

               <parameter key="ViscoBoundFactor" value="1" comment="Multiply viscosity value with boundary (default=1)" />

               <parameter key="DeltaSPH" value="0.1" comment="DeltaSPH value, 0.1 is the typical value, with 0 disabled (default=0)" />

               <parameter key="Shifting" value="0" comment="Shifting mode 0:None, 1:Ignore bound, 2:Ignore fixed, 3:Full (default=0)" />

               <parameter key="ShiftCoef" value="-6" 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="MPCoefficient" value="1" comment="Colagrossi and Landrini multiphase coefficient." />         

               <parameter key="RigidAlgorithm" value="1" comment="Rigid Algorithm 1:SPH, 2:DEM (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="DtIni" value="0.0000001" comment="Initial time step (default=h/speedsound)" units_comment="seconds" />

               <parameter key="DtMin" value="0.00000001" comment="Minimum time step (default=coefdtmin*h/speedsound)" units_comment="seconds" />

               <parameter key="#DtFixed" value="DtFixed.dat" comment="Dt values are loaded from file (default=disabled)" />

               <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="120" comment="Time of simulation" units_comment="seconds" />

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

               <parameter key="IncZ" value="2" comment="Increase of Z+" units_comment="decimal" />

               <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="0" comment="Minimum rhop valid (default=700)" units_comment="kg/m^3" />

               <parameter key="RhopOutMax" value="1300" comment="Maximum rhop valid (default=1300)" units_comment="kg/m^3" />         




  • And I post again missing pictures:

    t = 0 s

    t = 0.001 s

    t = 1.2 s

  • I hope someone with more experience with LiquidGas code can help find out if a bug is occuring or not - thanks for going into so much detail by the way.

    Hope it gets figured out, commenting to push this thread a bit higher up.

    Kind regards

  • Thanks for the idea, Asalih3D.

    I reallly don't understand how a subject like this doesn`t attracts more attention.

    Correct behaviour in simple workbenches like this is essential before proceeding with more complex analysis.

    This is my case. I am working in the development of a new kind of water column Wave Energy Converter. I have basic numeric models of it but a more detailed analysis is stalled. Behaviour of low compressibility is essential for me as low changes in air pressure are the core of the nehaviour in my system.

    Now I am refining my basic numeric model and I couldn't get advance in this subject.

    If I find more information about this subject I wil post it.

    Kind regardas

  • Dear @agavino

    I would suggest you to send an email to with more information about what your problem and implementation and we will try to help you or to forward your email to the member of the team that is now working in that direction


  • Hi, Alex.

    Thanks for the advice and sorry to have been offline so much.

    I already have sent e-mail. Wait for news.

    If I make some advance, I will post.


  • @agavino did you get news on this matter??

  • edited December 2019

    i could see you showed the equation:

    p = (cs0^2rhop0/gamma) ((rhop/rhop0)^gamma-1)

    but this is written in the guide for single phase. for multiphase, they say:

    "The multi-phase model uses a modified version of Tait’s equation of state [Batchelor, 1974] (the same equation used in single-phase DualSPHysics) for incompressible and inviscid fluids:


    i see your case is running only air, so maybe you should try it with the single phase code instead of two-phase code, as we run them separately. and when using two-phase code, im not sure neither where X (constant background pressure) comes from...

    anyway, i dont know why you got those results... if you got any answers, please share with us

  • im confused about that the EOS in two-phase code. according to the guide, X means backgound pressure, but i cannot see this term in the code.

    And, there is no any setting options in XML or argv. because X should be modifiable (by user). So, i suspect that it is not exist in dualsphysics.

    or,,, may be i just miss that.

  • @JOJO it does exist. You just need to add it in the source file where the pressure is defined, compile and run. In single phase, for GPU it is defined in

  • that is exactly what i do, it comfirms my suspect.

    But it is still werid that there is no X term in EOS ( or cpp,i mean it can not says it exists when one have to add it in code and compile again), because i think if the X term exist, it must be controled by some parameters or command, as a developer, one cannt let every user to find some terms or functions in code.Anywhere, i add a parameter to control it in xml.


  • Yeah, I agree. I don't think this is the most time efficient way. How did you add a parameter in .xml?

  • Very simple. Just like the other parameter ,such as mpcoe,h. Define a value in JSPH.cpp/.h and load it in xml(configconstant in JSphGpuSingle), then dont forget to updata it to Gpu in constantdataup in JSphGpu.cpp, then we can use it like ohter

  • i understood this case shouldnt be run in single-phase because weakly compressible approach would not fit gas behaviour properly and its compressibility.

    this adjust of X pressure, could you guys show me prints of how you´ve added it? @JOJO @circles

  • @agavino why did you use MPCoefficient=1?

    <parameter key="MPCoefficient" value="1" comment="Colagrossi and Landrini multiphase coefficient." />  

Sign In or Register to comment.