microstrain field

Asked by ytang

Hi all,
 I am trying to get the microstrain field. I did as the document said[1]. I want to use PyRunner command to export a series of .vtk files. It seems that I can just get one VTK file. is there someone can help me to check the ### get the strain field part ### in the MWE.

from yade import pack, plot
from yade import export
from yade import utils
compFricDegree=45.0
finalFricDegree=19.5
rate=-15
damp=0.7
target_strain = 0.15
stabilityThreshold=0.01
young=4e8
mn,mx=Vector3(0,0,0),Vector3(4e-3,8e-3,4e-3)
O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
psdSizes=[0.12e-3,0.14e-3,0.16e-3,0.18e-3,0.19e-3,0.20e-3,0.22e-3,0.24e-3,0.29e-3,0.38e-3]
psdCumm=[0,11.1,22.2,33.3,44.4,55.5,66.6,77.7,88.8,100]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=500,distributeMass=1,seed=1)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
triax=TriaxialStressController(
 maxMultiplier=1.0+2e5/young,
 finalMaxMultiplier=1.0+2e4/young,
 thickness=0,
 stressMask=7,
 internalCompaction=False
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 triax,
 newton,
]
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
 O.run(3000,True)
 unb=unbalancedForce()
 print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
 if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.01:
  break
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True
setContactFriction(radians(finalFricDegree))
re11=-triax.strain[0]
re22=-triax.strain[1]
re33=-triax.strain[2]
triax.stressMask = 5
newton=NewtonIntegrator(damping=0.0)
triax.goal2=rate
triax.strainRate2=15
triax.strainRate1=triax.strainRate3=1000.0
def stop_loading():
 if -triax.strain[1]-re22> target_strain:
  O.pause()
###################### get the strain field #####################################################
TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.setState(0)
#step = 0
#while -triax.strain[1]-re22< target_strain:
 #step = step+500
 #O.run(400,True)
 #TW.setState(1)
 #TW.defToVtk("strain_"+str(step)+".vtk")
O.run(100,True)
TW.setState(1)
def strain_export():
 #O.run(100,True)
 TW.setState(1)
 TW.defToVtk("strain.vtk")
#################################################################################################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 triax,
 newton,
 PyRunner(command='stop_loading()',iterPeriod=10000),
 PyRunner(command='strain_export()',iterPeriod=100),
]
O.run(50000,True)

I found that there is a question is almost the same as my question[2], they use a while function to solve this problem. But I want to know how can we solve this problem by using pyrunner function. thanks.

[1]https://github.com/yade/trunk/blob/master/examples/tesselationwrapper/tesselationWrapper.py
[2]https://answers.launchpad.net/yade/+question/675267

best,
yong

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
Last query:
Last reply:
Revision history for this message
ytang (ytang116) said :
#1

besides, when I run the code, the terminal pop up something I can not fully understand. the messages are as follows:

Triangulated Grains : 506
sym_grad_u_total_g (wrong averaged strain):
0.0425446 0.00274478 0.0282376
0.00274478 -0.211716 -0.0313596
0.0282376 -0.0313596 -0.0301356

Total volume = 1.79312e-09, grad_u =
0.0425446 -0.0116201 0.0640599
0.0171097 -0.211716 -0.0539018
-0.00758461 -0.00881733 -0.0301356

sym_grad_u (true average strain):
0.0425446 0.00274478 0.0282376
0.00274478 -0.211716 -0.0313596
0.0282376 -0.0313596 -0.0301356

Macro strain :
0 0 0
0 0 0
0 0 0
#################

what are the wrong average strain and true average strain?

thanks,

yong

Revision history for this message
Jan Stránský (honzik) said :
#2

Hello,

> It seems that I can just get one VTK file.
> TW.defToVtk("strain.vtk")

you specified only one file name, so you get only one file.
Have a look at [2] for an idea how to save multiple files. Will need some modification for PyRunner, something like:
###
class TWExporter:
   def __init__(self,base):
      self.base = base
      self.i = 0
   def export(self):
      TW.setState(1)
      TW.defToVtk("strain_{:03d}.vtk".format(self.i))
      self.i += 1
twExporter = TWExporter("strain")

O.engines = [
...
PyRunner(...,command="twExporter.export()"),
]
###

cheers
Jan

Revision history for this message
ytang (ytang116) said :
#3

Hi Jan,

So, from your answer, it means I need to create a class to contain the step (in your answer is self.i).

By the way, can you tell me what are the wrong average strain and true average strain printed out from the screen? what the relationship between these two values and the VTK.file that we exported from the export function?

thanks,
Yong

Revision history for this message
Jan Stránský (honzik) said :
#4

> So, from your answer, it means I need to create a class to contain the step (in your answer is self.i).

it is one option (there are other ways, so you do not need)

> can you tell me what are ...

no, sorry, I don't use and don't know these features..

Jan

Revision history for this message
Best Bruno Chareyre (bruno-chareyre) said :
#5

>what are the wrong average strain and true average strain printed out

Hi,
The true average strain is the proper volume average of local strain.
The wrong one is a weighted average where the weights do not follow from a proper partition of space.
Yet as you probably realize they are nearly identical.
Bruno

Revision history for this message
ytang (ytang116) said :
#6

Thanks Bruno Chareyre, that solved my question.

Revision history for this message
ytang (ytang116) said :
#7

Hi Jan,

thanks,

best,
yong

Revision history for this message
ytang (ytang116) said :
#8

Hi Bruno,

thank you!

Yong