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


Sign In or Register to comment.