ahout the pressure built-up on loose seabed

Hello all,

I'm trying to simulate the dynamic response of seabed under wave action by using DEM-PFV model.
But I have a problem when I try to numerical the loose seabed under wave action, the pore-pressure can't appear built-up and liquefaction occurred on seabed (refer to this in paper [1]). According to the paper[1] said, I create a loose seabed (flow.porosity=0.43), and increase the friction angle of sphere(frictionAngle=30). And I find that the fluid shear viscosity is 100 pa.s on the papare[1], but even set it , there is not built-up of pore-pressure occurring.
So I want to ask how to simulate this phenomenon ?
Here is my code:
***********************************************
                PREPARE SEABED
***********************************************
from yade import pack
from yade import *
import numpy as np
num_spheres=5000# number of spheres
young=1.5e7
compFricDegree = 20 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(2,1.75,0.75) # corners of the initial packing
pi=3.1415926

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2640,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalFricDegree),density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.1,num_spheres,periodic=True) #"seed" make the "random" generation always the same
spheres = [utils.sphere(s[0],s[1], material='spheres') for s in sp]
O.bodies.append(spheres)

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(gravity=(0,-9.81,0))

from yade import qt

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

## ______________ FLOW section _________________
# Activate flow engine and set boundary conditions
print "1"
triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
  O.run(100, True)
  unb=unbalancedForce()
  if unb<0.03 and abs(-10000-triax.meanStress)/10000<0.03:
    break
triax.stressMask=2
triax.goal1=triax.goal3=0
print '2'
triax.internalCompaction=False
triax.wall_bottom_activated=False
#load
triax.goal2=-11000; O.run(100,1)
print '3'
#unload
triax.goal2=-10000; O.run(100,1)
print "4"

O.save('compact-2,1.75,0.75-oedometer-30.yade')

***********************************************************
                     SIMULATE
***********************************************************

from yade import pack
from yade import *
from yade import qt
import numpy as np
O.timingEnabled=True
O.load('compact-2,1.75,0.75-oedometer-30.yade')
newton=NewtonIntegrator(gravity=(0,-9.81,0))

from yade import qt

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(label="flow"),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 #triax,
 newton
]
print "3"
print "n=",flow.porosity
flow.isActivated=1
flow.bndCondIsPressure=[0,0,0,1,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.meshUpdateInterval=10
flow.viscosity=0.1
flow.permeabilityFactor=1
flow.waveAction=1
flow.useSolver=3
flow.sineAverage=0
flow.slipBoundary=1
flow.tolerance=1e-06
flow.relax=1.9
flow.viscousShear=True
print "4"
t_zero=O.time
def timeDependentPressure():
       t = O.time-t_zero
       A =300
       T =1
       flow.sineMagnitude=A*sin(2*pi*t/T)
       flow.updateBCs()

O.engines=O.engines+[PyRunner(iterPeriod=20,command='timeDependentPressure()')]
##make nice animations:
O.engines=O.engines+[PyRunner(iterPeriod=100,command='flow.saveVtk()')]

O.saveTmp()
O.timingEnabled=1
from yade import timing

timing.stats()
*******************************************************************
Thank you for any suggestions.

Deng

[1]Catalano E. A pore-scale coupled hydromechanical model for biphasic granular media. An application to granular sediment hydrodynamics[J].