micro strain-paraview

Asked by Alireza Sadeghi on 2019-11-19

Hello all,

I want to simulate a sample in triaxial test. The sample has particles with three different sizes. The problem is that when I have only mono-size paricles, I can see the micro strain in paraview, but when I have three different size particles, the paraview does not show me the correct result. I was wondering if you could check my code.
Thank you very much.

Best Regards,

Alireza

####### code######

from yade import utils, plot
from yade import pack, qt
from datetime import datetime

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 0.31, fragile=False, label='aggregate-48'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1532.2, poisson=0.3, frictionAngle= 0.31, fragile=False, label='aggregate-814'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=2e11,
density=1523.6e2, poisson=0.3, frictionAngle= 0.28, fragile=False, label='wall'))

#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(2e-1,2e-1,2e-1)],thickness=0.0003,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t','t','t']

mats=['aggregate-48','aggregate-814','aggregate-814']

coke=((1.875e-3,500),(0.9475e-3,1367),(0.50125e-3,10618))

color=((0,0,1),(0,1,0),(1,1,0))

tolerance=[(2e-4),(1e-4),(1e-5)]

for i in range(len(nums)):

     nums[i]=pack.SpherePack()
     nums[i].makeCloud((0,0,0),(2e-1,2e-1,2e-1),rMean=coke[i][0],rRelFuzz=tolerance[i],num=coke[i][1])
     O.bodies.append([utils.sphere(c,r,material=mats[i],color=color[i]) for c,r in nums[i]])

#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
 [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
 [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
# VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
# qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
 NewtonIntegrator(damping=0.4,gravity=[0,0,-9.81])
]

yade.qt.Controller(), yade.qt.View()

O.dt=5e-8
#=================================================================
#================= APPLYING CONFINING PRESSURE ================
#=================================================================
sigmaIso=-1e5
triax=TriaxialStressController(

 maxMultiplier=1.000,
 finalMaxMultiplier=1.000,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True

while 1:
  O.run(100, 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 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb
  if triax.porosity<0.38 and abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.01:

    print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
    break
    print "### Isotropic state saved ###"

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('RVE-sizeDis-solid-Isoe5-Isopart.yade')

#====================================================================
#===================== DEVIATORIC LOADING =======================
#====================================================================

triax.stressMask = 3
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-1

O.saveTmp()

#====================================================================
#==================== Micro Strain ===========================
#====================================================================

TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
step=0
while 1:
    step +=1
    O.run(100,True)
    TW.setState(1)
    TW.defToVtk("strain_"+str(step)+".vtk")

O.engines=O.engines+[PyRunner(iterPeriod=20000,command='finishIt()',label='calmEngine')]

def finishIt():
 if abs(-triax.strain[2]) > 0.5:
# makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
  print O.iter
  print 'time is ', O.time
  print 'porosity:',triax.porosity
  print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2]
  O.save('RVE-sizeDis-solid-Isoe5-Devpart.yade')
  O.pause()

#=================================================================
#======================= data collecting =========================
#=================================================================

if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=100,command='history()',label='recorder')]+O.engines[5:7]

else:
  O.engines[4]=PyRunner(iterPeriod=100,command='history()',label='recorder')

def history():
   plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], ev=(triax.strain[0]+triax.strain[1]+triax.strain[2]), sxx=-triax.stress(triax.wall_right_id)[0], syy=-triax.stress(triax.wall_top_id)[1], szz=-triax.stress(triax.wall_front_id)[2], mass=sum(b.state.mass for b in O.bodies if isinstance(b.shape,Sphere)), q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, por=triax.porosity, i=O.iter),

plot.saveDataTxt('RVE-sizeDis-solid-Isoe5-Devpart.txt')

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-11-22
Last reply:
2019-12-08
Robert Caulk (rcaulk) said : #1

> the paraview does not show me the correct result

would you please describe the "result" you are encountering? Are there errors? If so what are the errors?

Alireza Sadeghi (asadeghime) said : #2

Dear Robert,

Thank you very much for your reply. When I open the micro strain.vtk files, the paraview give me these errors:

vtkVolumeTextureMapper3D was deprecated for VTK 7.0 and will be removed in a future version. vtkOpenGLVolumeTextureMapper3D was deprecated for V7X 7.0 will be removed in a future version.

and after it there is nothing in it. Does it mean that the problem is from my Paraview?I am using ParaView 5.4.1.
Thank you very much.

Best Regards,

Alireza

Robert Caulk (rcaulk) said : #3

I am happy to help you but it ran for 5 minutes trying to reach a porosity (before I quit the process, maybe it would've lasted for an hour? a day?). Is this long process of reaching a desired process a contributing factor to your issue? I tend to think not, can you please make a minimal working example (MWE) [1] that we can quickly and efficiently run to recreate your error?

[1]https://www.yade-dem.org/wiki/Howtoask

Robert Caulk (rcaulk) said : #4

desired porosity*

Alireza Sadeghi (asadeghime) said : #5

Dear Robert,

Thank you very much. I reduced the number of particles and now it takes time around 1 minute.
Sorry for the long code. In other cases, I don't have similar problem.
Thanks a lot.

Best Regards,

Alireza

############code#############

from yade import utils, plot
from yade import pack, qt
from datetime import datetime

qtr=qt.Renderer()
qtr.bgColor=(1,1,1)

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 0.31, fragile=False, label='aggregate-48'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1532.2, poisson=0.3, frictionAngle= 0.31, fragile=False, label='aggregate-814'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=2e11,
density=1523.6e2, poisson=0.3, frictionAngle= 0.28, fragile=False, label='wall'))

#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(2e-1,2e-1,2e-1)],thickness=0.0003,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t','t','t']

mats=['aggregate-48','aggregate-814','aggregate-814']

coke=((1.875e-3,50),(0.9475e-3,136),(0.50125e-3,106))

color=((0,0,1),(0,1,0),(1,1,0))

tolerance=[(2e-4),(1e-4),(1e-5)]

for i in range(len(nums)):

     nums[i]=pack.SpherePack()
     nums[i].makeCloud((0,0,0),(1e-1,1e-1,1e-1),rMean=coke[i][0],rRelFuzz=tolerance[i],num=coke[i][1])
     O.bodies.append([utils.sphere(c,r,material=mats[i],color=color[i]) for c,r in nums[i]])

#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
 [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
 [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
# VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
# qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
 NewtonIntegrator(damping=0.4,gravity=[0,0,-9.81])
]

yade.qt.Controller(), yade.qt.View()

O.dt=5e-7
#=================================================================
#================= APPLYING CONFINING PRESSURE ================
#=================================================================
sigmaIso=-1e5
triax=TriaxialStressController(

 maxMultiplier=1.000,
 finalMaxMultiplier=1.000,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True

while 1:
  O.run(100, 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 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb
  if triax.porosity<0.48 and abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.01:

    print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
    break
    print "### Isotropic state saved ###"

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('RVE-sizeDis-solid-Isoe5-Isopart.yade')

#====================================================================
#===================== DEVIATORIC LOADING =======================
#====================================================================

triax.stressMask = 3
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-1

O.saveTmp()

#====================================================================
#==================== Micro Strain ===========================
#====================================================================

TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
step=0
while 1:
    step +=1
    O.run(100,True)
    TW.setState(1)
    TW.defToVtk("strain_"+str(step)+".vtk")

O.engines=O.engines+[PyRunner(iterPeriod=20000,command='finishIt()',label='calmEngine')]

def finishIt():
 if abs(-triax.strain[2]) > 0.5:
# makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
  print O.iter
  print 'time is ', O.time
  print 'porosity:',triax.porosity
  print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2]
  O.save('RVE-sizeDis-solid-Isoe5-Devpart.yade')
  O.pause()

#=================================================================
#======================= data collecting =========================
#=================================================================

if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=100,command='history()',label='recorder')]+O.engines[5:7]

else:
  O.engines[4]=PyRunner(iterPeriod=100,command='history()',label='recorder')

def history():
   plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], ev=(triax.strain[0]+triax.strain[1]+triax.strain[2]), sxx=-triax.stress(triax.wall_right_id)[0], syy=-triax.stress(triax.wall_top_id)[1], szz=-triax.stress(triax.wall_front_id)[2], mass=sum(b.state.mass for b in O.bodies if isinstance(b.shape,Sphere)), q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, por=triax.porosity, i=O.iter),

plot.saveDataTxt('RVE-sizeDis-solid-Isoe5-Devpart.txt')

Robert Caulk (rcaulk) said : #6

After running your script (still not an MWE), the files open without error in paraview and show strain. Try downloading and reinstalling the latest paraview version.

Alireza Sadeghi (asadeghime) said : #7

Hello Dear Robert,

Thank you very much for your help. Actually the recent code (less particles) works great, and I can open the micro strain in paraview without any problem. But when I increased the number of particles as it was in the question (coke=((1.875e-3,500),(0.9475e-3,1367),(0.50125e-3,10618))), it does not work. My problem is "why the code does not work when I increase the number of particles?"
I tried it several times. Without any changes, when I increase the number of particles, I can't reach to micro strain.
Thank you very much again.

Best Regards,

Alireza

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

Hi, what's the meaning of «does not work» please?
Bruno

Alireza Sadeghi (asadeghime) said : #9

Hello Bruno,

When I increase the number of particles, as I said in the previeus massage, the strain tensor is not correct. In the code, I applied compaction strian rate (negetive) in z direction in the deviatoric part. Therefore, It should have negetive strain in the z direction. But now I have:

Triangulated Grains : 12491
sym_grad_u_total_g (wrong averaged strain):
0.0271423 -0.00121828 -0.00169358
-0.00121828 -0.00108962 -5.17863e-05
-0.00169358 -5.17863e-05 0.00136745

Total volume = 0.018195, grad_u =
0.0271423 6.67822e-06 -4.62223e-06
-0.00244324 -0.00108962 0.000223831
-0.00338253 -0.000327403 0.00136745

sym_grad_u (true average strain):
0.0271423 -0.00121828 -0.00169358
-0.00121828 -0.00108962 -5.17863e-05
-0.00169358 -5.17863e-05 0.00136745

Macro strain :
0 0 0
0 0 0
0 0 0

In addition, When I opened the micro strain in Paraview, there is no specific shape. However, when I open view in the controller, the shape of sample is correct. I don't know what should I do. All the conditions are the same, just I increased the number of particles.
Thank you very much.

best Regards,

Alireza

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

Can you reproduce the problem by prescribing motion to all bodies, like
vel=something*pos? Else it's unclear if the problem is in post processing
or in the simulation in itself.
Bruno

Le ven. 22 nov. 2019 00:09, Alireza Sadeghi <
<email address hidden>> a écrit :

> Question #686024 on Yade changed:
> https://answers.launchpad.net/yade/+question/686024
>
> Status: Answered => Open
>
> Alireza Sadeghi is still having a problem:
> Hello Bruno,
>
> When I increase the number of particles, as I said in the previeus
> massage, the strain tensor is not correct. In the code, I applied
> compaction strian rate (negetive) in z direction in the deviatoric part.
> Therefore, It should have negetive strain in the z direction. But now I
> have:
>
> Triangulated Grains : 12491
> sym_grad_u_total_g (wrong averaged strain):
> 0.0271423 -0.00121828 -0.00169358
> -0.00121828 -0.00108962 -5.17863e-05
> -0.00169358 -5.17863e-05 0.00136745
>
> Total volume = 0.018195, grad_u =
> 0.0271423 6.67822e-06 -4.62223e-06
> -0.00244324 -0.00108962 0.000223831
> -0.00338253 -0.000327403 0.00136745
>
> sym_grad_u (true average strain):
> 0.0271423 -0.00121828 -0.00169358
> -0.00121828 -0.00108962 -5.17863e-05
> -0.00169358 -5.17863e-05 0.00136745
>
> Macro strain :
> 0 0 0
> 0 0 0
> 0 0 0
>
> In addition, When I opened the micro strain in Paraview, there is no
> specific shape. However, when I open view in the controller, the shape of
> sample is correct. I don't know what should I do. All the conditions are
> the same, just I increased the number of particles.
> Thank you very much.
>
> best Regards,
>
> Alireza
>
> --
> 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
>
>
>

Alireza Sadeghi (asadeghime) said : #11

Dear Bruno,

Thank you very much for your reply.
I am not sure that the problem is in post processing. Because in low particles sample, the strain on the boundary is occurding to prescribed deformation, but when I increased the number of particles, the strain on the boundary does not follow the prescribed deformation. That is my problem.
What do you mean about prescribed motion to all bodies?Do you mean that I set a motion to each particle?What can it show?
thank you very much.

Best Regards,

Alireza

Launchpad Janitor (janitor) said : #12

This question was expired because it remained in the 'Open' state without activity for the last 15 days.