Yade

Contact laws and Input stiffness directly rather than young parameter

Asked by ehsan benabbas

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
ehsan benabbas
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hello,

> I want to directly input Kn and Kt ... instead of using the "young" parameter.

See [1,2,3] discussing similar topic

> To do so, as the "FrictMat" has not Kn and Kt, I tried to define material and contact law based on "ViscElMat".

you can use FrictMat, see [1], #3 and #8

> but, the result is different than the "FrictMat" case while the specimen is the same. ... I am wondering why that happens?

well, obviously, because you are using different model with different governing equations..

cheers
Jan

[1] https://answers.launchpad.net/yade/+question/679638
[2] https://answers.launchpad.net/yade/+question/298001
[3] https://answers.launchpad.net/yade/+question/228494

Revision history for this message
ehsan benabbas (ehsanben) said : 		#2

Hello

Thank you, Jan. As it's been mentioned in [1], this method leads to a SO LONG (more than 3 weeks) running time as I faced before. The problem of the long-run is the "young" parameter. When I calculate the "young" based on Kn, "young" will be a relatively high number, however, If I use "ViscElMat", the code with 20000 particles will be run less than 15 minutes. If I want to have the same run-time with "FrictMat" I need to forget about my desire Kn and use the value of "young" that Bruno already has used in his based code of triaxial test. I was wondering setting cn = cs = 0 in "ViscElMat" gives me the same result as "FrictMat" as the density of specimen has not changed (because of the specimen which is the same as when I used "FrictMat") and I should get the same behavior. For example, if I get the peak behavior (dense sample) with "FrictMat" maybe I get the same with "ViscElMat", which this did not happen.

[1] https://answers.launchpad.net/yade/+question/679638

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

> I was wondering setting cn = cs = 0 in "ViscElMat" gives me the same result as "FrictMat" as the density of specimen has not changed (because of the specimen which is the same as when I used "FrictMat") and I should get the same behavior.

why you "should get the same behavior", if you use different material, different Ip2 and different Law2?
And if you think / are sure that with certain parameters it leads to the same physics, why to change it?

> more than 3 weeks ... less than 15 minutes

I would check your stop conditions and which phase is "freezing" the simulation..

cheers
Jan

Revision history for this message
ehsan benabbas (ehsanben) said : 		#4

> why you "should get the same behavior", if you use a different material, different Ip2 and different Law2?
I get your point. However, it is not a different material (for me at least). When we are considering the same specimen, and this specimen is a real subject in our minds, we should get the same mechanical behavior apart from the scientific approach we look at the specimen. In fact, the material, the specimen, and its properties, all are the same, so, just changing the contact law should not lead to a different mechanical behavior (new stress-strain curve).

> And if you think/are sure that with certain parameters it leads to the same physics, why to change it?
With the FricMat, I get the dens behavior (with peak) which makes sense as the porosity is 0.4 and grains are in the sphere shape and fully rounded. The problem with this contact law is the "young" parameter which makes the code to be run so long as for 3 weeks. However, with the ViscElMat I get only the hardening behavior which refers to a loose specimen (with cn=cs=0) but the code will be run only in 15 minutes.
Again, for the FricMat, when I define the "young" based on Kn, it gets a value which makes the run-time so long.

Revision history for this message
ehsan benabbas (ehsanben) said : 		#5

> I would check your stop conditions and which phase is "freezing" the simulation.
My stop/finishing criterion is the target strain which is 0.25 (it's reasonable). With using FricMat and value for "young" in terms of Kn, The code will take so long in the deviatoric part and cannot reach the stop/finishing criterion easily.

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

> However, it is not a different material (for me at least).

well, different material objectively IS different material..

> just changing the contact law should not lead to a different mechanical behavior (new stress-strain curve).

sorry, but this is total nonsense. Actually, the main reason why to use a different contact law is to get different mechanical behavior..

> The problem with this contact law is the "young" parameter which makes the code to be run so long as for 3 weeks.
> However, with the ViscElMat I get only the hardening behavior which refers to a loose specimen (with cn=cs=0) but the code will be run only in 15 minutes.

you use the same young? then what are the differences in both approaches? Answering it should significantly help.

> Again, for the FricMat, when I define the "young" based on Kn, it gets a value which makes the run-time so long.

yes, that's one (not the only one) reason why this approach is discouraged (as discussed in the other referenced questions)

cheers
Jan

Revision history for this message
ehsan benabbas (ehsanben) said : 		#7

Thank you, Jan.

> you use the same young? then what are the differences in both approaches? Answering it should significantly help.
Actually I don't use "young" with the ViscElMat as I mentioned in the code. Let me copy that here:

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

I define these parameters and that is what I want as the input parameters (Kn, Ks, Density, Friction).

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

> Actually I don't use "young" with the ViscElMat
> However, it is not a different material (for me at least)

see?

> that is what I want as the input parameters (Kn, Ks, Density, Friction).

Do you know what is going on with ViscElMat instead of FrictMat? Me not, I have never used ViscElMat.
Maybe Kn is just a different name of a parameter but with the same meaning as young? You should be familiar with equations (Ip2, Law2) before further coparisons etc. (IMO).

cheers
Jan

Revision history for this message
ehsan benabbas (ehsanben) said : 		#9

I did something else, I think it's getting me the answer without adding a considerable time to the run-time. I defined a function as follows:

def stiff():
    for cont in O.interactions:
        cont.phys.kn = Whatever I want
        cont.phys.ks = Whatever I want

Then I defined a new PyRunner withe the frequency = 1 in "O.engines" as follows:

PyRunner(iterPeriod=1,command='stiff()',label='macro_recorder'),

Now I let the code to be run to see the result. I will talk about the result here for future reference. The only thing is I need to add a small code to check the value of Kn and Ks at the different parts of my code. What is your suggestion?

Revision history for this message
ehsan benabbas (ehsanben) said : 		#10

... new PyRunner with* the freque ...

PyRunner(iterPeriod=1,command='stiff()',label='stiff'),

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#11

The self-provided answer is correct in my view. It is an inelegant solution
to an inelegant problem. No offense to Ehsan, facing inelegant problems
happens to all of us.

For future refeference, better avoid such things.

Bruno

Le jeu. 13 févr. 2020 19:05, ehsan benabbas <
<email address hidden>> a écrit :

> Question #688681 on Yade changed:
> https://answers.launchpad.net/yade/+question/688681
>
> ehsan benabbas gave more information on the question:
> ... new PyRunner with* the freque ...
>
> PyRunner(iterPeriod=1,command='stiff()',label='stiff'),
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Revision history for this message
ehsan benabbas (ehsanben) said : 		#12

Bruno, thank you for your opinion, however, mine is totally different in the view of the research-side. The problem is the way Yade has been coded is not consistent with the most number of papers in micromechanics. To study micromechanics and the effect of contacts, researchers tend to have control over the stiffness values rather than the young parameter which is oversimplified. in this case more errors will be produced rather than the use of stiffness. In addition to that, researchers should be able to compare their works and results in order to make a conclusion. Otherwise, it is non-sense to only close the eyes and study without having a clear goal. This is why most of the users, particularly in micromechanics and geomechanics, like to use "Kn" and "Ks" as the input variables for contact law. in addition to that, such a long run-time is a serious obstacle in the way of researchers in this field. Therefore, paying attention to such issues of yade comes into the topic and this leads to using tricks instead of modifying the original code (as many similar cases have happened in FEM).

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#13

Could you please provide a blamelist of papers where, in 3D, stiffness is
defined in force/length instead of force/area?
Bruno

Le jeu. 13 févr. 2020 23:44, ehsan benabbas <
<email address hidden>> a écrit :

> Question #688681 on Yade changed:
> https://answers.launchpad.net/yade/+question/688681
>
> Status: Answered => Open
>
> ehsan benabbas is still having a problem:
> Bruno, thank you for your opinion, however, mine is totally different in
> the view of the research-side. The problem is the way Yade has been
> coded is not consistent with the most number of papers in
> micromechanics. To study micromechanics and the effect of contacts,
> researchers tend to have control over the stiffness values rather than
> the young parameter which is oversimplified. in this case more errors
> will be produced rather than the use of stiffness. In addition to that,
> researchers should be able to compare their works and results in order
> to make a conclusion. Otherwise, it is non-sense to only close the eyes
> and study without having a clear goal. This is why most of the users,
> particularly in micromechanics and geomechanics, like to use "Kn" and
> "Ks" as the input variables for contact law. in addition to that, such a
> long run-time is a serious obstacle in the way of researchers in this
> field. Therefore, paying attention to such issues of yade comes into the
> topic and this leads to using tricks instead of modifying the original
> code (as many similar cases have happened in FEM).
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

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

> for cont in O.interactions:
> cont.phys.kn = Whatever I want
> cont.phys.ks = Whatever I want

note that this leads to the same result as the solution proposed in [1] with constant E*r.
Furthermore, this is a worse solution, because you miss newly created interactions, which are evaluated within InteractionLoop with Law2 **before** this code is executed.

cheers
Jan

Revision history for this message
ehsan benabbas (ehsanben) said : 		#15

Yes, you are right Jan.
Just to mention, based on the formulation for Kn, we have Kn=E*r (if we assume E and r are constant for all particles)
While, based on the timestep calculation provided in the DEM formulation of the YADE website, in the section "Estimation of Δtcr
by wave propagation speed", this formula is Kn=E*D because π′ is equal to 2 in the classical DEM model (Ip2_FrictMat_FrictMat_FrictPhys) as implemented in Yade.

My conclusion is that as the particle mass, young, and size changes the timestep, it's not avoidable to have such a long runtime and the best solution is multiprocessing using as many cores as we can provide. Therefore, we can use E=Kn/D and need to think of a high-end PC or a supercomputer.

Thank you for your comments,

Bests,
Ehsan