Using flow engine got the wrong pressure field

Hello,

Thank you for kindly reading our posting guidelines [1]. Take care yo notice that we don't accept external links here.

As for your issue, it sounds like a stability problem, try decreasing the time step and turning O.dynDt=False.

[1]https://www.yade-dem.org/wiki/Howtoask

Cheers,

Robert

Hi Robert,

Thank you for answering my question. I'm sorry that I didn't read the posting guidelines carefully. I'll follow the guidelines in the future.

As for the issue, I tried decreasing the time step and turning
O.dynDt=False.
The result is almost the same. I noticed that the pressure field can be calculated by the following formula:

GP=Ev+Qq+Qp

where G is the conductivity matrix, P is the column vector containing all values of pressure, and E is the matrix defining the rates of volume change. Qq and Qp are flux vectors reflecting boundary conditions, respectively source terms (imposed fluxes) in Qq and imposed pressures in Qp.

I checked the initial velocity of my sample, it appears that all particles have a high velocity, thus in the second iteration gives a 'unreasonable' pressure field. I tried to slow down the particles, and the result is more and more 'reasonable'. Therefore, I think the way I generate the sample is improper. The code that I use to generate the sample is as followed:

####################################################
####################################################

from yade import pack

SoilMat=CohFrictMat(young=1e7,poisson=0.3,density=2650,frictionAngle=0.7,alphaKr=50,alphaKtw=50,momentRotationLaw=True)
O.materials.append((SoilMat))
O.materials.append(CohFrictMat(young=1e7,poisson=0.5,frictionAngle=0,density=0,label='walls',momentRotationLaw=True))

O.bodies.append(geom.facetBox((0,0,0.8),(.25,.25,0.8),wallMask=31,material='walls'))

sp=pack.SpherePack()
sp.makeCloud((-.25,-.25,0),(.25,.25,1.6),rMean=.0175,rRelFuzz=.3)
sp.toSimulation(material=SoilMat)

#Engine of the simulation
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 #introduced as a dead engine for the moment
 NewtonIntegrator(gravity=(0,0,-9.81)),
 #Run the specific controlled motion.
 PyRunner(command='start()',realPeriod=1,label='checker'),
]

O.dt=.5*PWaveTimeStep()

def start():
 newton.damping=0.3
 print('\n##### Simulation begins #####')
 print('Stage: Gravitity desposition.')
 print('Time step is %ss.' % O.dt)
 checker.command='checkUnbalanced()'

def checkUnbalanced():
 print('%s:UnbalancedForce=%s' % (O.iter,unbalancedForce()))
 if max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])>0.7: return
 if unbalancedForce()>0.01: return
 Initial_height=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
 print('The initial height of the simulation is %s' % Initial_height)
 print('Stage: Compression.')
 checker.command='Compression()'

def Compression():
 O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
 global plate
 plate=O.bodies[-1]
 plate.state.vel=(0,0,-.2)
 print('Wall successfully created! Initial veloctity is %s, initial position is %s.' % (plate.state.vel[2],plate.state.pos[2]))
 global Volume
 Volume=0
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   Volume=Volume+b.shape.radius*b.shape.radius*b.shape.radius*4*np.pi/3
 print('The volume of the pack is %s' % Volume)
 SoilMat.frictionAngle=0
 checker.command='Compression_S2()'

def Compression_S2():
 void_ratio=plate.state.pos[2]*.5*.5/Volume-1
 print('%s: Void ratio=%s' % (O.iter,void_ratio))
 if void_ratio>0.6: return
 SoilMat.frictionAngle=0.6
 plate.state.vel=(0,0,.2)
 print('Wait for the particles to settle down after compression.')
 checker.command='Density_Checker()'

def Density_Checker():
 void_ratio=plate.state.pos[2]*.5*.5/Volume-1
 print('%s: Void ratio=%s, Wall Position=%s' % (O.iter,void_ratio,plate.state.pos[2]))
 if abs(O.forces.f(plate.id)[2])>0: return
 if unbalancedForce()>0.01: return
 print(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]))
 checker.command='Record()'

def Record():
 print('The sample is generated successfully!')
 O.pause()

####################################################
####################################################

I think the way I generate the sample is not stable. And also, to spend less time, the material is different between when I generate the sample and when I do the seepage test.

I'm now trying to use the radius expansion method to generate the sample and use the same material in both test.
I hope this time it will work so that I can post my result as the end of this question.

Thanks Robert, it's indeed a stability problem.

I tried the radius expansion method, I think it's a fast way to generate stable sample. Anyway, both method (the old one I used before and the radius expansion method) can get the reasonable result in the end. It's just a matter of time.

I think the answer to my question is to find a way to generate particles that can reach force balance quickly.

Also, in order to get a stable pressure field, reducing the fluid viscosity helps a lot.

Creation of the specimen should have nothing to do with the stability of the coupled problem.

Your script is really unconventional and in some cases it is likely breaking yade. You are using start() as a PyRunner, which somehow continuously appends walls each time it is called. You should use oedometer.py [1] as a guide and rewrite your script.

[1]https://gitlab.com/yade-dev/trunk/blob/master/examples/FluidCouplingPFV/oedometer.py

Thanks, Robert.

I've rewrite my script. The particles generated by my old script appears to get a 'high' velocity in the first few iterations. And the particle velocity affects the change of pressure field. Therefore, I think careful creation of the specimen has certain effect on the change of pressure field. But the effect is little and it's not the key point. What I said not long ago was wrong.

Now I think the key point of the reason why the pressure field changes rapidly is the high value of the fluid viscosity. I change the viscosity from 10 to 0.00298 and found out that difference is obvious. The pressure field doesn't change that much.

>get a 'high' velocity in the first few iterations.

If this is true, then you just need to ensure that your unbalanced force is lower. This has nothing to do with flow engine.

>Now I think the key point of the reason why the pressure field changes rapidly is
>The pressure field doesn't change that much.

It changes rapidly or is doesnt change that much?

By the way, changing viscosity is the same as changing permeability in the code i.e. permeability and viscosity both contribute to intercell conductivity and therefore they have the same effect on stability and pressure field.

>If this is true, then you just need to ensure that your unbalanced force is lower. This has nothing to do with flow engine.

Yes, I agree with you.

>It changes rapidly or is doesnt change that much?

I found out that:

If the viscosity is high, the pressure field will change dramatically due to the external disturbance. Like my case, my unbalanced force is not low enough, so the first iteration gives pressure varies from 0-100 while the pressure of some pore is about -10000 in the second.

If the viscosity is low, the pressure field won't change much between the the first iteration and the second iteration. I think it's because I increase the value of the elements in the conductivity matrix G, therefore influence of the velocity of the particles is reduced.

GP=Ev+Qq+Qp

>changing viscosity is the same as changing permeability in the code

Thanks for the information. I didn't know this before.