yade-users team mailing list archive

[Rodrigo Borela <question270524@xxxxxxxxxxxxxxxxxxxxx>]

Question #270524 on Yade changed:
https://answers.launchpad.net/yade/+question/270524

Rodrigo Borela posted a new comment:
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 '========================================='

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.