2PFV - update Triangulation

Hi all,

I am investigating how the mesh evolves during the drainage using 2PFV.

I my original problem, I have a significant particle deformation/movement during the drying process, but tracking "getPoreThroatRadiusList()" in different instants I got the exactly same values. This also for the CellVolume, CellPorosity. It seems that the cell is constant even the pack is changing over time.

Any clue?

Here you have a simplified idea of my original script:

#!/usr/bin/python
# -*- encoding=utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/

from yade import pack
from yade import bodiesHandling
from yade import export
from yade import utils
from yade import ymport
import math
import numpy

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
    num_spheres=3000,# number of spheres
    compFricDegree = 1, # contact friction during the confining phase
    key='_triax_base_', # put you simulation's name here
    unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key

targetPorosity = 0.55 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0 # loading rate (strain rate)
damp=0.8 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
#2e4+70e4medio 1e4+70e4bom 1e4+60e4bom 3e4+90e4+w3,1,-1-the best
young=80e5 # contact stiffness200e4
young2=80e5
youngcoat=80e5
bondstr=1e3#2e7
bondstr2=1e3
bondstrcoat=1e6

## create materials for spheres and plates
mat=O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))
O.materials.append(JCFpmMat(type=1,young=20e7,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls'))
O.materials.append(JCFpmMat(type=1,young=youngcoat,poisson=0.3,frictionAngle=radians(1),density=1500,tensileStrength=bondstrcoat,cohesion=bondstrcoat,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstrcoat,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spherescoat'))

## create walls around the packing

mn,mx=Vector3(0,0,0),Vector3(0.0015,0.0015,0.0015)
mnbox,mxbox=Vector3(-0.0003,-0.0002,0.0003),Vector3(0.0025,0.0025,0.0025)#Vector3(0.002,0.00195,0.002)
walls=aabbWalls([mnbox,mxbox],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(minCorner=mnbox, maxCorner=mxbox,rMean=0.0001)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
#O.bodies.append(ymport.textExt("matrix_vtest6.txt", format='x_y_z_r', shift=Vector3(0,0.0001,0), scale=1.0,material='spheres',color=(0,1,1)))#0.00005
#O.bodies.append(ymport.textExt("coat_vtest6r756.txt", format='x_y_z_r', shift=Vector3(0,0,0), scale=1.0,material='spherescoat',color=(0,1,1)))
#O.bodies.append(ymport.textExt("coat_vtest6r56.txt", format='x_y_z_r', shift=Vector3(0,0,0), scale=1.0,material='spherescoat',color=(0,1,1)))

###########################
## DEFINING ENGINES ###
###########################

triax=TriaxialStressController(
    ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
    ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
    maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
    thickness = 0,
    ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
    ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
    ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
    ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
    ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask = 0,
    internalCompaction=False, # If true the confining pressure is generated by growing particles
    wall_front_activated=True,
    wall_back_activated=True,
    wall_top_activated=True,
    wall_bottom_activated=True,
    wall_left_activated=True,
    wall_right_activated=True,
    goal1=-200,
    goal2=-200,
    goal3=-200,
)

#newton=NewtonIntegrator(damping=damp)

O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1,label='Physicspheres')],#,xSectionWeibullShapeParameter=1.5, xSectionWeibullScaleParameter=1
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4),
 triax,
# VTKRecorder(iterPeriod=2000,recorders=['all','cracks'],fileName="/home/user/Área de Trabalho/Paper_Biopore/SandBox_Laser/Video/Videotest"),
# newton
# NewtonIntegrator(damping=damp,gravity=[0,-9.81,0]),
 newton

]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

## We will reach a prescribed value of porosity with the REFD algorithm
## (see http://dx.doi.org/10.2516/ogst/2012032 and
## http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

import sys #this is only for the flush() below

#poro=utils.voxelPorosityTriaxial(triax, 200, offset=0.0018)
while 1:
    O.run(50,True)
    # we decrease friction value and apply it to all the bodies and contacts
# compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print ("\r Friction: ",compFricDegree," porosity:",triax.porosity, triax.stress(3)[1],triax.stress(0)[0],triax.stress(4)[2],triax.strain[0],triax.strain[2],
    sys.stdout.flush())
    # while we run steps, triax will tend to grow particles as the packing
    # keeps shrinking as a consequence of decreasing friction. Consequently
    # porosity will decrease
    unb=unbalancedForce()
# triax.goal1=triax.goal2=triax.goal3=triax.meanStress*2
    print(unb)
    if triax.strain[0]<-0.10:
      triax.goal1=0
      triax.goal3=0
    if triax.porosity<targetPorosity:
# triax.goal1=triax.goal2=triax.goal3=0
      break

#O.save('compactedState'+key+'.yade.gz')
print ("### Compacted state saved ###")
print(triax.stress(3)[1])

##############################
### DEVIATORIC LOADING ###
##############################

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
#triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))
#O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2000,tensileStrength=bondstr2,cohesion=bondstr2,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr2,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))

#=======================================

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 2
triax.wall_bottom_activated=1
#now goal2 is the target strain rate
triax.goal1=rate#triax.stress(1)[0]
triax.goal3=rate#triax.stress(5)[2]
triax.goal2=triax.stress(3)[1]

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
#O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot
from yade import plot
O.run(10,True)

#strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L (strain)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

# a function saving variables
def history():
    plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            s11=-triax.stress(triax.wall_right_id)[0]-si0,
            s22=-triax.stress(triax.wall_top_id)[1]-si1,
            s33=-triax.stress(triax.wall_front_id)[2]-si2,
            e=math.exp(-triax.strain[1]-ei1)-1,
            pc=-unsat.bndCondValue[2],
            sw=unsat.getSaturation(isSideBoundaryIncluded=False),
            z1=O.bodies[3].state.pos[1],
            i=O.iter)

if 1:
  # include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

plot.plots={'pc':('sw',None,'e22')}
plot.plot()
print("test",math.exp(2))

#######################################################
## Drainage Test under oedometer conditions ###
#######################################################
##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
#meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondIsWaterReservoir=[0,0,1,0,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=False #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension = 0.0728

#start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt

ts=O.dt
pgstep= 40#45000000*ts #30Pa/s
print (pgstep)
pgmax= 10000#9316 #Pa
mi=0.0009 #Pa.sec
f1=open('cells.txt',"w")

for pg in arange(1.0e-8,pgmax,pgstep):
  print(pg)
  unsat.bndCondValue=[0,0,(-1.0)*pg,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()

  unsat.permeabilityMap=True

  dy=utils.aabbDim()
  q=max(unsat.averageVelocity())
  L=dy[1]#*(1+triax.strain[1]+ei1)
  P=(mi*q*L)/pg

# print(celsV1)
# print(unsat.getSaturation(isSideBoundaryIncluded=False),pg,-triax.strain[1],dy[1],q,L,P)

  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))

  unsat.meshUpdateInterval=20
  unsat.defTolerance=-1
  unsat.updateTriangulation=True
  unsat.breakControlledRemesh=True

# if pg==520.00000001:
  cels=unsat.nCells()
# celsV1=[0.0]*cels
  for ii in range(cels):
    celsV1=numpy.array(unsat.getCellCenter(ii))
    celsV2=numpy.array(unsat.getCellVelocity((celsV1[0],celsV1[1],celsV1[2])))
    f1.write(str(celsV2)+"\n")
  f1.close()

  while 1:
    O.run(100,True)
# unsat.breakControlledRemesh=True
# unsat.updateTriangulation=True
    unb=unbalancedForce()
    if unb<0.1:
      break