yade-users team mailing list archive

[Leonard <question688791@xxxxxxxxxxxxxxxxxxxxx>]

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

Leonard posted a new comment:
Hi Jan,
Thanks for your suggestion.
Here are the two scripts based on Yade 2018.02b:
#########
The following is the script1 which is for generating sample (takes around 30 seconds):
# -*- coding: utf-8 -*-
#*************************************************************************
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import timeit

##########  this script is modified from https://gitlab.com/yade-
dev/trunk/blob/master/examples/triax-tutorial/script-session1.py#L142

utils.readParamsFromTable(lowerR=10.0,upperR=10.0)
from yade.params import table

num_spheres=2000# number of spheres
targetPorosity = 0.5 #the porosity we want for the packing
compFricDegree = 30 # 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.01 # loading rate (strain rate)
damp=0.4 # damping coefficient
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)
young=5e7 # contact stiffness

confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

## create materials for spheres and plates
MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(CohFrictMat(isCohesive=False,young=young,poisson=0.5,frictionAngle=radians(30),
										 density=2650.0,normalCohesion=1e6, shearCohesion=1e6,label='sand'))
## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3333,num_spheres,True, 0.95,seed=1)

O.bodies.append([sphere(center,rad,material='sand') for center,rad in
sp])

Gl1_Sphere.quality=3

triax=TriaxialStressController(
	maxMultiplier=1.+6e5/young, # spheres growing factor (fast growth)
	finalMaxMultiplier=1.+2e4/young, # spheres growing factor (slow growth)
	thickness = 0,
	stressMask = 7,
	internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
		 Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
	newton,
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
	break


print "###    Isotropic completed     ###"

import sys
while triax.porosity>targetPorosity:
	# 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,
	sys.stdout.flush()

        O.run(500,1)

print "###  state 2  Reach target porosity completed      ###"

print 'top', -triax.stress(triax.wall_top_id)[1],
'right',-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2]

triax.internalCompaction=False
triax.stressMask = 7
triax.goal1=triax.goal2=triax.goal3=-confinement
triax.max_vel=0.001

O.save('sample.yade.gz')
#####################################################
The following is script2 which is for triaxial test, it illustrates the problem.
##
# -*- coding: utf-8 -*-
#*************************************************************************
from yade import pack, plot
import math
import numpy as np

O.load('sample.yade.gz')
Gl1_Sphere.quality=3
confinement=100e3

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
	thickness = 0,
	stressMask = 5,
	internalCompaction=False, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.4)


O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
		 Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
    PyRunner(command='stopIfDamaged()',iterPeriod=1000),
	newton,
]

Gl1_Sphere.stripes=0
yade.qt.Controller()
yade.qt.View()
setContactFriction(radians(30))
triax.stressMask = 5
triax.goal2=-0.015
triax.goal1=-confinement
triax.goal3=-confinement

from yade import plot

### a function saving variables
def history():
	plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
			Volumetric_strian=triax.strain[0]+triax.strain[1]+triax.strain[2],
			s11=-triax.stress(triax.wall_right_id)[0],
			s22=-triax.stress(triax.wall_top_id)[1],
			s33=-triax.stress(triax.wall_front_id)[2],
			dev=-triax.stress(triax.wall_top_id)[1]-confinement,
			i=O.iter)

def stopIfDamaged():
	if -triax.strain[1]>0.15:
		O.pause()
		plot.saveDataTxt('data')


if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+[PyRunner(iterPeriod=100,command='stopIfDamaged()',label='checkDamage')]+O.engines[5:7]

plot.plots={'e22':('s11','s22','s33'),'e22 ':('dev'),'e22  ':('Volumetric_strian')}
plot.labels={'e22':'$\epsilon_{a}$','s11':'$\sigma_{11}$','s22':'$\sigma_{22}$','s33':'$\sigma_{33}$'}
plot.plot()
########
Run script2 after script1, I have:
In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured: 
InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines.

You may switch the position of newton with PyRunner in the O.engines (in
script2) to see the difference. For me, the error is solved by this way.

Cheers,

Leonard

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.