yade-users team mailing list archive

[ehsan benabbas <question688681@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello everyone,

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

I developed the Triaxial test code based on [1] and working on that. I want to directly input Kn and Kt (like most of the papers) instead of using the "young" parameter. To do so, as the "FrictMat" has not Kn and Kt, I tried to define material and contact law based on "ViscElMat". To do so and as I do not want to model visco materials, I set cn = cs = 0 as follows:

O.materials.append(ViscElMat(kn=Kn,ks=Kt,cn=0.0,cs=0.0,density=particleDensity,frictionAngle=radians(compFricDegree),label='spheres'))
O.materials.append(ViscElMat(kn=Kn*100,ks=Kt*100,cn=0.0,cs=0.0,density=0,frictionAngle=0,label='walls'))

and for the contact law in Engine, I made the following changes:

newton=NewtonIntegrator(damping=damp)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
	[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
        [Law2_ScGeom_ViscElPhys_Basic()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
    PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
	TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key),
	newton
]

In this case, I get no error, the run speed is fast enough, but, the result is different than the "FrictMat" case while the specimen is the same. To be more specific, I get a dense behavior (Peak in the stress-strain curve) by using "FrictMat" however, when I use the "ViscElMat" I get a loose soil behavior (just gradually hardening) and never use "young" in my code. I am wondering why that happens? 

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

Thank you for your help.

Minimal version of my code:

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

import numpy as np
import time
import datetime, os
import shutil
start_time=time.time()
from datetime import datetime
import math
from yade import qt, export, utils
from yade import pack
import yade.timing;
O.timingEnabled=True

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

nRead=readParamsFromTable(
 num_spheres=20000,
 compFricDegree = 29,
 key='_triax_',
 unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.4
compFricDegree = table.compFricDegree
finalFricDegree = 30
IP=100 # iteration period to record data and stuff
ORN=1000 # O.Run Number of iterations
damp=0.2
thick=0.01
stabilityThreshold=0.01
PCPC=0.0001     # Precision of Confining Pressure Convergence
Kn=10e7
Kt=10e7
Ls=1  
mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls)
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
cvs=1-targetPorosity    # coefficient of solid volume ===>>> 1 - n = 1 - (Vv / Vt) = Vs / Vt   , cvs*volume=Vs
mean_rad = pow(0.24*cvs*volume/num_spheres,0.3333)
strainRate=-0.02 
target_strain=0.25
sigmaIso=-1e4 
particleDensity=2600

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

O.materials.append(ViscElMat(kn=Kn,ks=Kt,cn=0.0,cs=0.0,density=particleDensity,frictionAngle=radians(compFricDegree),label='spheres'))
O.materials.append(ViscElMat(kn=Kn*100,ks=Kt*100,cn=0.0,cs=0.0,density=0,frictionAngle=0,label='walls'))

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

walls=aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
clumps=False
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

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

triax=TriaxialStressController(
 maxMultiplier=1.005,
 finalMaxMultiplier=1.002,
 thickness = thick,
 stressMask = 7,
 internalCompaction=True,
)

##################################################
################# 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,
        t = O.time,     # virtual (yade) time --- time of simulation
        fab = utils.fabricTensor()[0])

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

newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
	[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
        [Law2_ScGeom_ViscElPhys_Basic()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
    PyRunner(iterPeriod=IP,command='history()',label='macro_recorder'),
	TriaxialStateRecorder(iterPeriod=IP,file='WallStresses'+table.key),
	newton
]

Gl1_Sphere.stripes=True
yade.qt.Controller(), yade.qt.View()

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

triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = sigmaIso

while 1:
    O.run(ORN,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
    ConfStressRatio=abs(sigmaIso-triax.meanStress)/abs(sigmaIso)
    if unb<stabilityThreshold and ConfStressRatio<PCPC:
        break

e22Check=-triax.strain[1]

########################################################
################# REACHING TARGET POROSITY ##################

import sys
while triax.porosity>targetPorosity:
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	sys.stdout.flush()
	O.run(ORN,True)

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

triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))

triax.wall_top_activated=True
triax.wall_bottom_activated=False
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True

triax.stressControl_1=True
triax.stressControl_2=False
triax.stressControl_3=True

triax.strainRate2=strainRate

triax.stressMask = 5
triax.goal1 = sigmaIso
triax.goal2 = strainRate
triax.goal3 = sigmaIso

newton.damping=0.1

while 1:
    O.run(ORN,True)
    unb=unbalancedForce()
    axialS=triax.stress(triax.wall_top_id)[1]
    eps2=triax.strain[1]
    eps2cal=(triax.height-Ls)/Ls
    if abs(eps2-e22Check)>=target_strain: #triax.sigma2
        break

##########################
########### END ###########

O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')
print('Real time of run:',(time.time() - start_time), 'seconds or',(time.time() - start_time)/60, 'minutes')
yade.timing.stats()

#######################################################
################# RECORD Macro DATA and Plot #################

from matplotlib import pyplot as plt

e22 = plot.data["e22"]
ev = plot.data["ev"]
s22 = plot.data["s22"]
s11 = plot.data["s11"]
s33 = plot.data["s33"]
i = plot.data["i"]
t = plot.data["t"]

fig1, ax1=plt.subplots(figsize=(15, 10))
fig2, ax2=plt.subplots(figsize=(15, 10))
fig3, ax3=plt.subplots(figsize=(15, 10))
fig4, ax4=plt.subplots(figsize=(15, 10))

ax1.plot(e22,ev,color='k',linewidth=1,label='Volumetric Strain')
ax1.legend()
ax1.set_xlabel('Axial Strain')
ax1.set_ylabel('Volumetric Strain')
ax1.grid(True)
fig1.tight_layout()
fig1.savefig('plot1.eps', format='eps')

ax2.plot(t,s22,color='k',linewidth=1,label='Axial Stress')
ax2.legend()
ax2.set_xlabel('Time')
ax2.set_ylabel('Axial Stress')
ax2.grid(True)
fig2.tight_layout()
fig2.savefig('plot2.eps', format='eps')

ax3.plot(e22,s22,color='k',linewidth=1,label='Stress-Strain')
ax3.legend()
ax3.set_xlabel('Axial Strain')
ax3.set_ylabel('Axial Stress')
ax3.grid(True)
fig3.tight_layout()
fig3.savefig('plot3.eps', format='eps')

ax4.plot(t,s11,color='b',linewidth=1,label='Confining Stress - 1')
ax4.plot(t,s33,color='r',linewidth=1,label='Confining Stress - 3')
ax4.legend()
ax4.set_xlabel('Time')
ax4.set_ylabel('Confining Pressure')
ax4.grid(True)
fig4.tight_layout()
fig4.savefig('plot4.eps', format='eps')

plt.show()

O.pause()


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