yade-users team mailing list archive

[ehsan benabbas <question687838@xxxxxxxxxxxxxxxxxxxxx>]

New question #687838 on Yade:
https://answers.launchpad.net/yade/+question/687838

Hi everyone,

I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74

I use the Triaxial code by Bruno Chareyre [1] to understand how it works.

[1] https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py

My goal is to code a Triaxial test with 20000 particles in different sizes (sand) under different strain and stress paths and get all inter-particle forces, branch vectors, contact normals, and fabric tensor. (with the properties defined in the code)

This is the code I am running. My questions are at the end of this message.

###################################################
############################################################################################################################
#########         	TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS         #########
############################################################################################################################

import datetime
aa = datetime.datetime.now()
print ('************** START **************')
import numpy as np
import math
from yade import pack, plot, qt, export, utils
from datetime import datetime

######################################################
#########         DEFINING VARIABLES         #########

print ('============ DEFINING VARIABLES ============')
nRead=readParamsFromTable(
	num_spheres=20000,
	compFricDegree = 30,
	key='_triax_base_',
	unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres		
key=table.key
targetPorosity = 0.35 			
compFricDegree = table.compFricDegree 	
finalFricDegree = 30 			
rate=-0.02 				
damp=0.2 				
stabilityThreshold=0.01 		
young=5e6 			
poisson=0.3
mn,mx=Vector3(0,0,0),Vector3(20,20,20)
sigmaIso=-50e3				

######################################################
#########         DEFINING MATERIALS         #########

print ('============ DEFINING MATERIALS ============')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))		
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionlesswalls'))				

####################################################
#########         DEFINING PACKING         #########

print ('============ DEFINING PACKING ============')
frictionlesswalls=aabbWalls([mn,mx],thickness=0,material='frictionlesswalls')				
wallIds=O.bodies.append(frictionlesswalls)									
sp=pack.SpherePack()								
clumps=False 											
if clumps:
    volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])						
    mean_rad = pow(0.09*volume/num_spheres,0.3333)							
    c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])	
    sp.makeClumpCloud(mn,mx,[c1],periodic=False)							
    sp.toSimulation(material='spheres')								
    O.bodies.updateClumpProperties()								
else:
    volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
    sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) 					
    O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

##########################################################
#########         DEFINING TRIAXIAL TEST         #########

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(			
	maxMultiplier=1.+2e4/young, 		
	finalMaxMultiplier=1.+2e3/young, 	
	thickness = 0,
	stressMask = 7,
	internalCompaction=True,		
)
newton=NewtonIntegrator(damping=damp)

####################################################
#########         DEFINING ENGINES         #########

print ('============ DEFINING ENGINES ============')
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),	
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
	newton
]
Gl1_Sphere.stripes=True											
if nRead==0: yade.qt.Controller(), yade.qt.View()
print ('Number of elements: ', len(O.bodies))
print ('Box Volume: ',  triax.boxVolume)
print ('Box Volume calculated: ', volume)

###############################################################
#########         APPLYING CONFINING PRESSURE         #########

print ('============ APPLYING CONFINING PRESSURE ============')		
triax.stressmask=7
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
while 1:
    O.run(1000, True)
    unb=unbalancedForce()	
    meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3	
    print ('unbalanced force:',unb,' mean stress engine: ',triax.meanStress,' mean stress (Calculated): ',meanS)
    print ('porosity=',triax.porosity)
    print ('void ratio=',triax.porosity/(1-triax.porosity))
    print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
    if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
        break
O.save('confinedPhase'+key+'.xml')
print ('################## Isotropic phase is finished and saved successfully ##################')
e22Check=triax.strain[1]

######################################################
#########         DEVIATORIC LOADING         #########

print ('============ APPLYING DEVIATORIC LOADING ============')	
triax.internalCompaction=False	
setContactFriction(radians(finalFricDegree))	
triax.stressMask = 5						
triax.goal1=sigmaIso
triax.goal2=rate				
triax.goal3=sigmaIso	
newton.damping=0.1		
unb=unbalancedForce()	
axialS=triax.stress(triax.wall_top_id)[1]
print ('step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-sigmaIso)
print ('axial deformation (%)', (triax.strain[1]-e22Check)*100)
print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
O.save('final.xml')
print ('################## Deviatoric phase is finished and saved successfully ##################')

######################################################
#########         DEFINING FUNCTIONS         #########

print ('============ DEFINING FUNCTIONS ============')
def history():											
	plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
			ev=-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],
			i=O.iter)

########################################################
#########         RECORD AND PLOT DATA         #########

print ('============ RECORD AND PLOT DATA ============')
O.run(5000,True)
plot.plots={'e22':('s11','s22','s33',None,'ev')}						
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$'}
plot.plot()											
plot.saveDataTxt('results'+key)								
plot.saveGnuplot('plotScript'+key)
print ('************** END **************')
bb = datetime.datetime.now()
cc = bb - aa
print( int(cc.total_seconds()))
print(elapsed_time)
###################################################

Question:
1- When the simulation is finished? When I run the code, It takes some minutes to show the ************** END ************** message on the screen, while the Play button has not been chosed in the Yade window yet. If I press the Play button, it takes some hours and then I get run errors.
2- In the code we have defined the strain rate for the deviatoric part, what's the criterion to finish the test? does it work with damping? I have read about it but it's not clear for me yet.
3- I want to get the time of running but I get the error for it. I have used datetime function for it as it is in the code.
4- I don't get the conventional graph of triaxial test. the graph window shows up but it is just a point.
5- If I get the simple triaxial test correctly, how can I get the fabric tensor, inter-particle forces, branch vectors, and contact normals as the raw data?

Thank you.

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