yade-users team mailing list archive

Thread
Date
[Question #687871]: Outputs of the Triaxial code by Bruno Chareyre

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: ehsan benabbas <question687871@xxxxxxxxxxxxxxxxxxxxx>
Date: Thu, 09 Jan 2020 01:11:54 -0000
Reply-to: question687871@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
New question #687871 on Yade:
https://answers.launchpad.net/yade/+question/687871

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)

Based on the code and its output which can be seen in the following, I have asked my questions at the end of this message.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

print ('************** START **************')
import numpy as np
import datetime
import time
import math
from yade import 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 FUNCTIONS         #########

print ('============ DEFINING FUNCTIONS ============')
from yade import plot
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,
			n = triax.porosity)

######################################################
#########         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)

from yade import pack								
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,		
)

####################################################
#########         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,
	PyRunner(iterPeriod=20,command='history()',label='recorder'),
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
	NewtonIntegrator(damping=damp,label="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 ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')    
    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))
    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 ##################')

########################################################
#########         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 **************')
O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now output:

ehsan@ehsan:~/Desktop$ /home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 Q.py
Welcome to Yade 2019-08-08.git-775ae74 
Using python version: 3.6.9 (default, Nov  7 2019, 10:44:02) 
[GCC 8.3.0]
TCP python prompt on localhost:9000, auth cookie `scdyua'
XMLRPC info provider on http://localhost:21000
Running script Q.py
************** START **************
============ DEFINING VARIABLES ============
============ DEFINING FUNCTIONS ============
============ DEFINING MATERIALS ============
============ DEFINING PACKING ============
============ DEFINING TRIAXIAL TEST ============
============ DEFINING ENGINES ============
The constructor with a shareWidget is deprecated, use the regular contructor instead.
Number of elements:  20006
Box Volume:  6.90392657805617e-310
Box Volume calculated:  8000.0
============ APPLYING CONFINING PRESSURE ============
~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~
unbalanced force: 0.007284146090760937  mean stress engine:  -1.3899805976256234  mean stress (Calculated):  -2.537923714531715
porosity= 0.8343868539188172
void ratio= 5.038168005756043
################## Isotropic phase is finished and saved successfully ##################
============ APPLYING DEVIATORIC LOADING ============
step= 1000 unbalanced force: 0.007284146090760937  sigma2:  -0.0 q= 50000.0
axial deformation (%) 0.0
~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~
################## Deviatoric phase is finished and saved successfully ##################
============ RECORD AND PLOT DATA ============
/home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/plot.py:444: MatplotlibDeprecationWarning: 
The 'verts' kwarg was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use 'marker' instead.
  scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,s=scatterSize,color=line.get_color(),**scatterMarkerKw)
************** END **************
Analysis has been taken for 270.053 seconds or 4.500883333333333 minutes
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

In [1]:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Questions:

1- Is there any problem with getting the following message in the output right after "DEFINING ENGINES"?
>> The constructor with a shareWidget is deprecated, use the regular contructor instead.

2- Why the volume box calculated by triax is near to 0? 
>> Box Volume:  6.90392657805617e-310

3- Why in the "APPLYING CONFINING PRESSURE" section the porosity is about 0.8 while I defined it as 0.35?
>> porosity= 0.8343868539188172

Thank you,
Ehsan

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