Energy of the body is very high, defying gravity and normal stress - cpm material

Asked by Akm on 2020-02-01

Hi all,

I am trying to simulate a shear test where I have two CPM materials with different properties - top_mat(lower strength) and bot_mat(higher strength). top_mat rests on bot_mat at an angle <45. Try to imagine a direct shear test where there is an asperity at the interface. There is only friction between the two materials-no cohesion.

There is a plate at the top which applies a load on top_mat in Constant Normal stiffness condition(CNS). top_mat has its boundaries free only in the Y direction where it is allowed to move vertically up or down. bot_mat can only move in the X direction(horizontal).

When bot_mat(bottom box) is moved using a plate in the horizontal direction(X), the top_mat should ''SLIDE'' along the asperity in the Y direction and then be subjected to damage when the shear strength threshold is reached for top_mat.

But even though I have a plate with a CNS load applied on top, during the shearing process, top_mat starts to ''FLY'' up instead of sliding and loses contact with bot_mat by pushing the CNS plate up, thereby defying both the CNS and the gravitational forces in the system. I am not sure why this is happening.

P.S: I had modified the source code to have zero cohesion ONLY AT THE INTERFACE: iscohesive=0, undamagedcohesion=0

Sorry that I am unable to think of a way to provide an MWE because I cut out that interface using a gts surface and this platform doesn't support the inclusion of external links. If it is alright, I can provide the code and the external links here. Kindly help me out as I have been stuck in this issue for many days. Help as soon as possible will be much appreciated. Many thanks in advance.

-Arun

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-02-27
Last reply:
2020-02-27

This question was reopened

  • 2020-02-25 by Akm
Akm (arunkumarmurali) said : #2

I will post the full script here. I will attach the drive links for the particles because I exported those particles using a gts surface. I did the sample preparation in a separate file using radial expansion. So the particle positions are here-kindly download those text files.

''https://drive.google.com/drive/folders/1M2rKs-vnR3uHC_ybv2vliYLhvO2g0RTm?usp=sharing''

The complete script is as follows:

from yade import pack, qt, export, ymport, plot

#Default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
 intfactor = 1.2, #Max value of 1.2 is preferred as it can be reverted back to 1 after a step.
     Shear_Vel=0.1,
 S_length=0.060,
 S_area=0.060*0.030,
 Sig0=300e3, #300kN/m2
 kn=600e3, #300kPa/mm

 plotperiod=500,
 topname='/home/akm/Documents/Yade/install/bin/Full-10k-Top.txt',
 botname='/home/akm/Documents/Yade/install/bin/Full-10k-Bot.txt',
 spname='/home/akm/Documents/Yade/install/bin/snapshots-trial/plot',
 vtkname='/home/akm/Documents/Yade/install/bin/vtk-trial/s',
)

from yade.params.table import *
F_init=Sig0*S_area #Initial Normal force
#MATERIAL SPECIFICATION

mat1=CpmMat(
 young = 32e9,
     density=2400,
 poisson = 0.30,
 frictionAngle = radians(35),
 epsCrackOnset = 5e-1,
 relDuctility = 1.1,
 sigmaT = 1.2e6,
)

mat2=CpmMat(
 young = 5e9,
     density=1940,
 poisson = 0.30,
 frictionAngle = radians(35),
 epsCrackOnset = 5e-2,
 relDuctility = 1.0001,
 sigmaT = 1.2e1,
)
O.materials.append(mat1)
O.materials.append(mat2)

#UPPER PARTICLE IMPORT:
upbox_particles=ymport.text(topname,material=mat2, color=(1,0,0))
O.bodies.append(upbox_particles)

#LOWER PARTICLE IMPORT:
lowbox_particles=ymport.text(botname,material=mat1, color=(0,1,0))
O.bodies.append(lowbox_particles)

###################
##GEOMETRY CREATION :
####################

O.materials.append(FrictMat(frictionAngle=0,density=0,label='walls'))

#Dimension of Packing
d1,d2=Vector3(.01,.01,.01),Vector3(.07,.0783,.04) # corners of the initial packing

xmin=d1[0]
xmax=d2[0]
ymin=d1[1]
ymax=d2[1]
zmin=d1[2]
zmax=d2[2]
ht_lower=0.03
ht_upper=0.03
asp_ht=0.0083

#Shear plate for shearing in the +X direction
P_Shear=utils.box((xmin,ymin+(ht_lower/2)-(ht_lower/20),zmin+(zmax-zmin)/2 ), (0,(ht_lower/2),(zmax-zmin)/2 ),fixed=True, wire=True, color=(0,0,1))
O.bodies.append(P_Shear)
P_Shear.state.blockedDOFs = "xyzXYZ"

#Surrounding setup geometry
U_box_1=geom.facetBox((xmin+(xmax-xmin)/2, ymin+ht_lower+asp_ht+(asp_ht+ht_upper)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht+ht_upper)/2,(zmax-zmin)/2 ), wallMask=2) # +X Boundary constraints
U_box_2=geom.facetBox((xmin+(xmax-xmin)/2, ymin+ht_lower+(asp_ht*2+ht_upper)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht*2+ht_upper)/2,(zmax-zmin)/2 ), wallMask=1) # -X Boundary constraints
L_box=geom.facetBox((xmin+(xmax-xmin)/2, ymin+(asp_ht+ht_lower)/2,zmin+(zmax-zmin)/2 ), ((xmax-xmin)/2,(asp_ht+ht_lower)/2,(zmax-zmin)/2 ), wallMask=6)
Cover_box=geom.facetBox(((xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2 ),((xmin+xmax+0.06)/2,(ymin+ymax+0.02)/2,(zmin+zmax+0.02)/2 ) ,wallMask=63)
O.bodies.append(U_box_1)
O.bodies.append(U_box_2)
O.bodies.append(L_box)
O.bodies.append(Cover_box)

idlist=[]
idlist.append(P_Shear.id)
for i in L_box:
    idlist.append(i.id) #ids for the translation motion of the box

for i in O.bodies:
 if isinstance(i.shape,Sphere):
  if i.mat==mat1:
   i.state.blockedDOFs = "yzXYZ" #blocking the mat1 movement to ensure that there is no moment at the end point

#Adding a plate at the top:
global plate
plate_pos=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if isinstance(b.shape,Sphere)])
plate=utils.box((xmin+(xmax-xmin)/2,plate_pos,zmin+(zmax-zmin)/2), ((xmax-xmin)/2 - xmin/60, 0, (zmax-zmin)/2),fixed=True, wire=False,color=(0,0,1))
O.bodies.append(plate)

print('Plate added')
plate.state.blockedDOFs = "xyzXYZ"

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intfactor,label='bo1s'),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intfactor,label='ig2s'),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)],
  [Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom_CpmPhys_Cpm()]
 ),
    NewtonIntegrator(gravity=(0,-9.81,0),damping=0.5),
    GlobalStiffnessTimeStepper(active=True,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8),
 PyRunner(iterPeriod=1,command='addPlotData()'),
]

O.step()
bo1s.aabbEnlargeFactor=1
ig2s.interactionDetectionFactor=1

def addPlotData():

 Sf = -O.forces.f(P_Shear.id)[0] #Negative sign since the force applied is in -X direction
 Hor_dspl = 1000*P_Shear.state.displ()[0] #Hor disp in mm
 Shear_stress = Sf/S_area
 Vert_dspl = 1000*plate.state.displ()[1]
 Force_Stif=kn*Vert_dspl*S_area #Force(N)=Stiffness(N/m2/mm)*Displacement(mm)*Area(m2)
 Force_Reqd=Force_Stif+F_init #Final force on plate = Force due to CNS + Intial force on the plate(SigmaO)
 Forceonplate=O.forces.f(plate.id)[1]
 if Forceonplate<Force_Reqd:
  plate.state.vel=(0,-2,0)
 elif Forceonplate>Force_Reqd:
  plate.state.vel=(0,1,0)
 Norm_stress_applied=(Forceonplate)/S_area
 Norm_stress_calculated=Force_Reqd/S_area

 yade.plot.addData(
  Hor_Disp = Hor_dspl,
  Shear_stress = Shear_stress,
  Ver_Disp = Vert_dspl,
  Norm_stress_applied = Norm_stress_applied,
  Norm_stress_calculated=Norm_stress_calculated
 )

# note the space in 'i ' so that it does not overwrite the 'i' entry
plot.plots={'Hor_Disp':('Shear_stress'),'Hor_Disp ':('Norm_stress_applied','Norm_stress_calculated'),}
plot.plot()

#EXTRA ENGINES:
O.engines=O.engines+[TranslationEngine(ids=idlist,translationAxis=[1,0,0], velocity=Shear_Vel)]
O.engines=O.engines+[SnapshotEngine(fileBase=spname,iterPeriod=plotperiod)]
#O.engines=O.engines+[VTKRecorder(fileName=vtkname,recorders=['cpm','all'],iterPeriod=500)]

yade.qt.View()
O.saveTmp()

#####Thanks in advance guys#####

Jan Stránský (honzik) said : #3

Hello,

what version of Yade you use?

> top_mat
> bot_mat

next time please make the question and code consistent

> There is only friction between the two materials-no cohesion.

not in the provided code

> But even though I have a plate with a CNS load applied on top, during the shearing process, top_mat starts to ''FLY'' up instead of sliding and loses contact with bot_mat by pushing the CNS plate up, thereby defying both the CNS and the gravitational forces in the system. I am not sure why this is happening.

I could not reproduce this (probably because there is cohesive connection between top and bottom)
But, according to
> top_mat has its boundaries free only in the Y direction where it is allowed to move vertically up or down
then your description seems logical to me.. Is
Shear_Vel=0.1m/s
realistic? Maybe a real specimen would also jump in reality moving horizontally 10 cm per second (just guessing what could be wrong).

> P.S: I had modified the source code to have zero cohesion ONLY AT THE INTERFACE: iscohesive=0, undamagedcohesion=0

please provide the code

cheers
Jan

Akm (arunkumarmurali) said : #4

Hi Jan,

Thanks for looking into the code.

>what version of Yade you use?

It is 1.2.0. For the sake of modifyng the source code, I had downloaded the Yade/trunk files - https://github.com/yade/trunk

 and then installed Yade using three folders as outlined in the installation manual:

https://www.yade-dem.org/doc/installation.html?highlight=sudo%20add%20apt%20repository%20ppa%20yade%20users%20external

So when I run it, I get it as Yade Unknown.

>next time please make the question and code consistent

Sorry, next time, will do.

>I could not reproduce this (probably because there is cohesive connection between top and bottom)

I have only changed three lines in the source code. It will be not good to paste the entire thing here. Let me put that source code in the drive as well so that you can download it and compile it.

''https://drive.google.com/drive/folders/1M2rKs-vnR3uHC_ybv2vliYLhvO2g0RTm?usp=sharing''

>realistic? Maybe a real specimen would also jump in reality moving horizontally 10 cm per second

We know that the physical nature of the experiment cannot be simulated here in the simulations. This is a 10000 particle sample. 0.1mps is a 125% strain rate which covers like 5% strain in 20 minutes. I tried to use a 12.5% strain rate as well, it covered 5% of the strain in 13 hours. But the sample behaved the same way like it did with higher strain rates. There is absolutely no difference.

Jan Stránský (honzik) said : #5

>>what version of Yade you use?
> It is 1.2.0

1.20.0?

> I have only changed three lines in the source code.

in this case, you could have post the three lines instead of complete file..

> We know that the physical nature of the experiment cannot be simulated here in the simulations.

then you cannot expect physical nature of the results..
you can try some simple hand calculations, maybe you just imposed such conditions that real gravity and plate pressure simply does not have enough effect..

cheers
Jan

Akm (arunkumarmurali) said : #6

Hi Jan,

It's 1.20.0.

Were you able to figure out any discrepancy in my code that led to the overall behaviour?

The conditions are valid and experiments have been performed under those conditions. I was looking for an approximate behaviour like that. Please let me know if you found any issues out here.

Many thanks,
Arun

Jérôme Duriez (jduriez) said : #7

Yade 1.20.0 is extremely old (before June 2016), see https://launchpad.net/yade/+download

Jan Stránský (honzik) said : #8

> It's 1.20.0.
> I have only changed three lines in the source code.

Comparing the provided file ConcretePM.cpp and the one from 1.20.0, there is much more than three lines changes..

> Were you able to figure out any discrepancy in my code that led to the overall behaviour?
> epsCrackOnset = 5e-1,
> epsCrackOnset = 5e-2

Where these values come from?
this epsCrackOnset would lead to non-realistically high strength, very roughly f=E*e=32e9*5e-1=16 GPa / 1,6 GPa. Then the material behaves as elastic and you will not get any damage..
The elastic sample cannot do anything else (damage is not possible and stiffness is too high to somehow deform dramatically) than to ''FLY'' up..

Jan

Akm (arunkumarmurali) said : #9

I think that might have been the problem. I had given the tensile strength as too high and that has created a lot of problems for me. How do I know which version of Yade I am using? I had downloaded the package from https://github.com/yade/trunk.
 So it shows the version as Yade unknown. Should I be downloading it from somewhere else?

Many thanks,
Arun

Jan Stránský (honzik) said : #10

> How do I know which version of Yade I am using?

start Yade, it should print something like "Welcome to Yade XYZ"

Jan

Jérôme Duriez (jduriez) said : #11

And forget about github:

See "The main [Yade] branches [= source code] are currently on GitLab (migrated from GitHub in January 2019)." on https://yade-dem.org/doc/

Akm (arunkumarmurali) said : #12

Sure many thanks.

Akm (arunkumarmurali) said : #13

Sorry to bring this problem up again. I thought the 'flying up' phenomenon was due to very high tensile strength. But even when I use strengths of 0.8Mpa, I still can see that the sample rebounds and flys up, thus losing contact with the bottom box. I could verify it in two ways :

1. Seeing the images in ParaView - the upper box was literally up in the air without any contact with the lower box.
2. Mathematically, the residual shear strength should be "Normal stress*tan(phi)". But the residual went to zero.

Even when I do the same code for just Direct shear where I have a plain interface, there is the same behaviour. The top sample is losing contact with the bottom sample at any cost. Can't we model the direct shear of two rigid bodies using Cpmmat?

I tried using JCFPm as well but could see the same results. Can someone please guide me into this?

Thanks,
Arun

Jan Stránský (honzik) said : #14

What happen if you increase the force on the plate 10x? 1000x? 1000000x?
cheers
Jan

Akm (arunkumarmurali) said : #15

Hi Jan,

I increased the normal load by 10 times, but still, the residual strength drops to zero. If I increase by 100x or something, the entire sample crushes under the load and there is no point in doing the simulation for that high normal stress. Please help me understand the cause and solution to this problem.

Thanks,
Arun

Jan Stránský (honzik) said : #16

what version of Yade are you actually using (see above discussion)?

> Please help me understand the cause and solution to this problem.

there might be millions of reasons..
E.g., if I got it correctly, you prescribe position of the loading plate according to displacement and reaction based on a stiffness value.
Are you sure the stiffness is correct? (just seeing mixing m and mm, even if it is correct, I suggest to use basic SI units or at least using one unit for one quantity)
The plate apart from stiffness has some mass which is omitted. This may influence the results?
...

Resume, it is just a model and maybe it just does not suit your needs, or maybe there are some assumptions that are not realistic, or maybe some values are wrong (stiffness, velocity, ...), or maybe the material has non-realistic properties, or ........

cheers
Jan

Akm (arunkumarmurali) said : #17

Hi,

I downloaded the source code from gitlab - it mentions the yade as Yade unknown while running.

I have checked everything that you said - but something is not right and I am unable to figure it out.

Chareyre (bruno-chareyre-9) said : #18

Please read installation chapter of documentation first, and if you have
installation issues report in another question with appropriate details of
what you did.

Also please avoid questions related to some source code you modified.
Especially if the changes are not disclosed.

Finally, please, make the scripts minimal and problem statement short if
possible.
That's to increase your chances of getting usefull feedback.

Bruno

Le jeu. 27 févr. 2020 10:57, Akm <email address hidden> a
écrit :

> Question #688400 on Yade changed:
> https://answers.launchpad.net/yade/+question/688400
>
> Status: Needs information => Open
>
> Akm gave more information on the question:
> Hi,
>
> I downloaded the source code from gitlab - it mentions the yade as Yade
> unknown while running.
>
> I have checked everything that you said - but something is not right and
> I am unable to figure it out.
>
> --
> 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
>
>
>

Can you help with this problem?

Provide an answer of your own, or ask Akm for more information if necessary.

To post a message you must log in.