Evolving material properties and timestep instability

Hello Rodrigo,
please provide a MWE [1], without it it is very difficult to guess what is
the reason of the problem
thanks
Jan

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

2015-08-19 6:11 GMT+02:00 Rodrigo Borela <
<email address hidden>>:

> New question #270524 on Yade:
> https://answers.launchpad.net/yade/+question/270524
>
> Dear all,
>
> I am working on the modelling of a problem using the JCFpm contact law on
> a very dense pack (porosity is about 11%), where I have three phenomena
> being imposed simultaneously:
> 1. particle shrinkage, applied using utils.growParticles(shrinkage factor
> of the order of .9999995)
> 2. stiffening of the contacts, applied using an iterative loop on O.bodies
> (b.mat.young*=stiffeningfactor) and an iterative loop on O.interactions (
> i.phys.kn*=stiffeningfactor and i.phys.ks*=stiffeningfactor)
> (stiffeningfactor of the order of 1.0000005)
> 3. increasing tensile strength, applied also using iterative loops on
> O.bodies and O.interactions following the same procedure as in (2).
> I'm using the GlobalStiffnessTimeStepper, where the iteration interval (of
> about 500-1000 iterations) is the same as the one for when the three
> phenomena are applied.
> The problem I'm experiencing is that the time step has been decreasing at
> very large rates despite the gentle stiffening+shrinking, for example from
> 1e-8 to 1e-10 in an iteration interval. At some point, the particles
> explode or the time step becomes 'nan'.
> The work I'm doing is reproducing results from a study that employed PFC,
> and all of my coefficients and ratios are calculated according to their
> work.
> Has anyone experienced something similar and would have a suggestion to
> handle this seemingly numerical instability?
>
> Thank you in advance for your time!
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Dear Jan,

thank you for your reply. Here follows the code. The text file containing my sample information is here http://paste.ubuntu.com/12129242/

O.reset()
import math, sys, os
from yade import utils, pack, plot, ymport, timing

O.timingEnabled=False

# Total experiment duration (iterations)
tau=1000000.0

#--------- Initial Material properties
# Friction angle
phi1=0
tanPhi1=0
# Tensile strength
sigmat1=3.729e6
# Young's modulus
young1=1.25e7

# Add material to the simulation
cohFric1=O.materials.append(JCFpmMat(cohesion=sigmat1,tensileStrength=sigmat1,frictionAngle=phi1,young=young1,density=2790,poisson=1,label='cohFric1'))

# -------------- Cylindrical Sample
# Center
centerx=0.0
centery=0.0
centerz=0.0
# Dimensions
segNumber=60
sampleDiameter=0.12612
h=0.01484
top=centerz+(h/2)
bottom=centerz-(h/2)
# Cylinder
O.bodies.append(geom.facetCylinder(center=(centerx,centery,centerz),radius=sampleDiameter/2,height=h,orientation=Quaternion((1,0,0),0),segmentsNumber=segNumber,wallMask=6))
# Particle pack
O.bodies.append(ymport.text('cylindrical-initialD18.75-24687part.txt'))

#=====================================================================
# ENGINES
#=====================================================================

#-------------- Iteration interval
nIter=int(tau/1000.0)

# Engines
O.engines=[
         ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key='test-rborelav',recordCracks=True)],
   ),
   PyRunner(command='EvolvingPropertiesAndShrinkage()',iterPeriod=nIter),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
   GlobalStiffnessTimeStepper(active=True,dead=False,timeStepUpdateInterval=nIter)
]

#=====================================================================
# EVOLVING PROPERTIES AND SHRINKAGE
#=====================================================================

#---------------- Evolution stages
 # Initial
wi=0.60
# Intermediate 1
w1=0.17
# Intermediate 2
w2=0.11
# Final
wf=0.045

# Shrinkage parameter
alpha=.198
# Stiffness parameter
B=-0.197
# Tensile strength parameter
D=-0.083

w=wi
def EvolvingPropertiesAndShrinkage():
      global nIter,tau,wi,w1,w2,wf, w, alpha
      global B, D
      # Current iteration
      t=O.iter
      # Water content
      w=((t/tau)*(wf-wi))+wi
      precedingW=(((t-nIter)/tau)*(wf-wi))+wi
      # Stiffness multiplier
      stiffnessMultiplier=(exp(B*w))/(exp(B*precedingW))
      # Tensile strength multiplier
      sigmatMultiplier=(exp(D*w))/(exp(D*precedingW))
      #------------- Evolving material properties
      # Material properties update
      for b in O.bodies:
            # Normal stiffness
            b.mat.young*=stiffnessMultiplier
            # Tensile strength
            b.mat.tensileStrength*=sigmatMultiplier
            # Cohesive part of shear strength
            b.mat.cohesion*=sigmatMultiplier
      # Interactions update
      for i in O.interactions:
            # Normal stiffness
            i.phys.kn*=stiffnessMultiplier
            # Shear stiffness
            i.phys.ks*=stiffnessMultiplier
            # Tensile Strength
            i.phys.FnMax*=sigmatMultiplier
            # Cohesive part of shear strength
            i.phys.FsMax*=sigmatMultiplier
      #------------- Shrinkage
      shrinkRatio=(exp(-alpha*(t/tau)))/(exp(-alpha*((t-nIter)/tau)))
      for b in O.bodies:
            utils.growParticles(shrinkRatio)
      #------------- Information
      print 'The current time step is ', O.dt
      print 'The current water content is ',w
      print 'The average number of contacts is ',utils.avgNumInteractions
      print 'The shrink ratio is ',shrinkRatio
      print 'The stiffness multiplier is ',stiffnessMultiplier
      print 'The tensile strength multiplier is ',sigmatMultiplier
      print '========================================='

Sorry! I forgot a pair of brackets in front of "utils.avgNumInteractions". I corrected it on this one. Many thanks for the patience!

O.reset()
import math, sys, os
from yade import utils, pack, plot, ymport, timing

O.timingEnabled=False

# Total experiment duration (iterations)
tau=1000000.0

#--------- Initial Material properties
# Friction angle
phi1=0
tanPhi1=0
# Tensile strength
sigmat1=3.729e6
# Young's modulus
young1=1.25e7

# Add material to the simulation
cohFric1=O.materials.append(JCFpmMat(cohesion=sigmat1,tensileStrength=sigmat1,frictionAngle=phi1,young=young1,density=2790,poisson=1,label='cohFric1'))

# -------------- Cylindrical Sample
# Center
centerx=0.0
centery=0.0
centerz=0.0
# Dimensions
segNumber=60
sampleDiameter=0.12612
h=0.01484
top=centerz+(h/2)
bottom=centerz-(h/2)
# Cylinder
O.bodies.append(geom.facetCylinder(center=(centerx,centery,centerz),radius=sampleDiameter/2,height=h,orientation=Quaternion((1,0,0),0),segmentsNumber=segNumber,wallMask=6))
# Particle pack
O.bodies.append(ymport.text('cylindrical-initialD18.75-24687part.txt'))

#=====================================================================
# ENGINES
#=====================================================================

#-------------- Iteration interval
nIter=int(tau/1000.0)

# Engines
O.engines=[
         ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key='test-rborelav',recordCracks=True)],
   ),
   PyRunner(command='EvolvingPropertiesAndShrinkage()',iterPeriod=nIter),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
   GlobalStiffnessTimeStepper(active=True,dead=False,timeStepUpdateInterval=nIter)
]

#=====================================================================
# EVOLVING PROPERTIES AND SHRINKAGE
#=====================================================================

#---------------- Evolution stages
 # Initial
wi=0.60
# Intermediate 1
w1=0.17
# Intermediate 2
w2=0.11
# Final
wf=0.045

# Shrinkage parameter
alpha=.198
# Stiffness parameter
B=-0.197
# Tensile strength parameter
D=-0.083

w=wi
def EvolvingPropertiesAndShrinkage():
      global nIter,tau,wi,w1,w2,wf, w, alpha
      global B, D
      # Current iteration
      t=O.iter
      # Water content
      w=((t/tau)*(wf-wi))+wi
      precedingW=(((t-nIter)/tau)*(wf-wi))+wi
      # Stiffness multiplier
      stiffnessMultiplier=(exp(B*w))/(exp(B*precedingW))
      # Tensile strength multiplier
      sigmatMultiplier=(exp(D*w))/(exp(D*precedingW))
      #------------- Evolving material properties
      # Material properties update
      for b in O.bodies:
            # Normal stiffness
            b.mat.young*=stiffnessMultiplier
            # Tensile strength
            b.mat.tensileStrength*=sigmatMultiplier
            # Cohesive part of shear strength
            b.mat.cohesion*=sigmatMultiplier
      # Interactions update
      for i in O.interactions:
            # Normal stiffness
            i.phys.kn*=stiffnessMultiplier
            # Shear stiffness
            i.phys.ks*=stiffnessMultiplier
            # Tensile Strength
            i.phys.FnMax*=sigmatMultiplier
            # Cohesive part of shear strength
            i.phys.FsMax*=sigmatMultiplier
      #------------- Shrinkage
      shrinkRatio=(exp(-alpha*(t/tau)))/(exp(-alpha*((t-nIter)/tau)))
      for b in O.bodies:
            utils.growParticles(shrinkRatio)
      #------------- Information
      print 'The current time step is ', O.dt
      print 'The current water content is ',w
      print 'The average number of contacts is ',utils.avgNumInteractions()
      print 'The shrink ratio is ',shrinkRatio
      print 'The stiffness multiplier is ',stiffnessMultiplier
      print 'The tensile strength multiplier is ',sigmatMultiplier
      print '========================================='

Thanks Jan Stránský, that solved my question.

Thank you very much, Jan. I changed that and ran the simulation from beginning to end and it worked perfectly :)