In WCSPH, why does the static water pressure continue to rise after pressure balance?

[hrbyrb]

I made an interesting calculation example these days. I used WCSPHSheme to simulate the process of water standing. During this process, I observed that the pressure of water fluctuated greatly, and after almost reaching the pressure balance, the pressure of water would continue to rise for a long time.

The video is shown below:

dam_break_3d_WCSPH_hydro.mp4

My code:

import numpy as np
from pysph.base.kernels import WendlandQuintic
from pysph.examples._db_geometry import DamBreak3DGeometry
from pysph.solver.application import Application
from pysph.sph.integrator import EPECIntegrator
from pysph.sph.scheme import WCSPHScheme

# parameter to change the resolution
dim = 3
dt = 1e-4
tf = 5.0
dx = 0.05
nboundary_layers = 1
hdx = 1.3
ro = 1000.0
h0 = dx * hdx
gamma = 7.0
alpha = 0.25
beta = 0.0
c0 = 10.0 * np.sqrt(2.0 * 9.81 * 0.55)


class DamBreak3D(Application):
    def add_user_options(self, group):
        group.add_argument(
            '--dx', action='store', type=float, dest='dx',  default=dx,
            help='Particle spacing.'
        )
        group.add_argument(
            '--hdx', action='store', type=float, dest='hdx',  default=hdx,
            help='Specify the hdx factor where h = hdx * dx.'
        )

    def consume_user_options(self):
        dx = self.options.dx
        self.dx = dx
        self.hdx = self.options.hdx
        self.geom = DamBreak3DGeometry(
            container_height=1.0, container_width=1.0, container_length=1.2+dx,
            fluid_column_height=0.55, fluid_column_width=1.0, fluid_column_length=1.2,
            obstacle_center_x=2.5, obstacle_center_y=0,
            obstacle_length=0.16, obstacle_height=0.161, obstacle_width=0.4,
            dx=dx, nboundary_layers=nboundary_layers, hdx=self.hdx, rho0=ro,
            with_obstacle=False
        )
        self.co = 10.0 * self.geom.get_max_speed(g=9.81)

    def create_scheme(self):
        s = WCSPHScheme(
            ['fluid'], ['boundary'], dim=dim, rho0=ro, c0=c0,
            h0=h0, hdx=hdx, gz=-9.81, alpha=alpha, beta=beta, gamma=gamma,
            hg_correction=True, tensile_correction=False, delta_sph=False
        )
        return s

    def configure_scheme(self):
        s = self.scheme
        hdx = self.hdx
        kernel = WendlandQuintic(dim=dim)
        h0 = self.dx * hdx
        s.configure(h0=h0, hdx=hdx)
        #dt = 0.25*h0/(1.1 * self.co)
        s.setup_properties(self.particles)
        s.configure_solver(
            kernel=kernel, integrator_cls=EPECIntegrator, tf=tf, dt=dt,
            adaptive_timestep=False, n_damp=100 ,pfreq=100,
            output_at_times=[0.4, 0.6, 1.0]
        )

    def create_particles(self):
        return self.geom.create_particles()

    def customize_output(self):
        self._mayavi_config('''
        viewer.scalar = 'p'
        b = particle_arrays['boundary']
        b.plot.actor.mapper.scalar_visibility = False
        b.plot.actor.property.opacity = 0.1
        ''')

if __name__ == '__main__':
    app = DamBreak3D(output_dir="F:/sph/outputfile/dambreak_test_WCSPHScheme_1g")
    app.run()
    app.post_process(app.info_filename)