← Back to team overview

yade-users team mailing list archive

Re: [Question #270524]: Evolving material properties and timestep instability

 

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

Jan Stránský proposed the following answer:
Hi Rodrigo,

thanks for the script

      for b in O.bodies:
>             # Normal stiffness
>             b.mat.young*=stiffnessMultiplier


b.mat is the same material instance for all bodies, so yout code would mean
that you multiply mat.young *= pow(stiffnessMiltiplier, numberOfBodies),
making mat.young very high number. Changing material properties does not
influence existing contacts (so the problem might be visible much later),
but when new contact is created, it could be the problem (as you mentioned,
big reduction of time step by GlobalStiffnessTimestepper)

Instead, try something like:

mats = set(b.mat for b in O.bodies)
for mat in mats:
   mat.young *= stiffnessMultiplier
   ...

please try this approach and let us know if it helped or not
cheers
Jan



2015-08-19 20:12 GMT+02:00 Rodrigo Borela <
question270524@xxxxxxxxxxxxxxxxxxxxx>:

> Question #270524 on Yade changed:
> https://answers.launchpad.net/yade/+question/270524
>
> Rodrigo Borela posted a new comment:
> 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 '========================================='
>
> --
> 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp
>

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