# Contact Particle + Nmesh and Nmesh + triangulated.box

17 views

Hi

'm simulating railroad ballast, and for that, the lab tests applied the cyclic load on the dormant, which transfer to the ballast. The process is to insert the particules and compress them, so the dormant is inserted after the particules and, in the end, the box that apply the cyclic load, through a harmonic oscillation, is inserted. I need to know if there is something I can do to make contact between this materials.

1. Between the dormant (Nmesh import) and the triangulated box
2. Between the dorment and the particules, is 1 can't be done

First I tried 2, but even when I tried to reduce dtSafety from 0,8 to 0,5, when the dormant touch the particles, it explod.

Is the contact of these materials implemented or I need to uso something else? Or am I doing something wrong?

I'm open to suggestion.

Regards
Alan

from woo.core import *
from woo.dem import *
from minieigen import *
from woo.fem import *
import woo.triangulated
from math import pi
import woo.models, woo

woo.master.usesApi=10102

# ioneMask=2 >> Particle groups which have bits in loneMask in common (i.e. (A.mask & B.mask & loneMask)!=0) will not have contacts between themselves

# contact model plus material for particles
model=woo.models.ContactModelSelector(name='linear',damping=.4,mats=[FrictMat(young=5e5,ktDivKn=.2,tanPhi=.4)])
# same material for walls as for particles, but without friction
S.lab.wallMat=model.mats[0].deepcopy(tanPhi=0.)
# intial box size
xmaxbox = 2.
ymaxbox = 2.
xdor = .62
ydor = .22
zdor = .15
box=AlignedBox3((0,0,0),(xmaxbox,ymaxbox,4))
lastro=AlignedBox3((0,0,1.5),(2,2,3.5))
slastro=AlignedBox3((0,0,0),(2,2,2))
# use something else than spheres here, if you like
gen2=PsdSphereGenerator(psdPts=[(.2,0),(.25,1)])
gen=PsdSphereGenerator(psdPts=[(.1,0),(.2,1)])
ctrlStep=100
m0=woo.utils.defaultMaterial()
m1=m0.deepcopy()
m0.young*=5
m0.density=1e4 # make it faster

# wall box without the top (will be added later)

S.engines=DemField.minimalEngines(model=model)+[
# generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
# check what's happening periodically
PyRunner(ctrlStep,'checkProgress(S)')]+[IntraForce([In2_Tet4_ElastMat(),In2_Facet(),In2_Membrane_FrictMat(bending=True)])]

S.lab.stage=0
# avoid warning when we change the value
S.lab._setWritable('stage')
# things that we want to look at during the simualtion
S.plot.plots={'dz':('Fz0','Fz1')}

S.saveTmp()

# callback function from PyRunner
def checkProgress(S):
# if the top wall is there, save its position and acting force
if S.lab.stage in (1,2):
# graphic adjustment z0 x Fz0
elif S.lab.stage in (4,5):
# graphic adjustment z1 x Fz1
if S.lab.stage==0 and S.lab.inletL.dead and woo.utils.unbalancedForce(S)<0.05:
S.lab.stage+=1
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
z0=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
# place Wall just above the top
S.lab.topWall=Wall.make(z0,axis=2,sense=-1,mat=S.lab.wallMat)
# print z0 before compression
print('z0' , z0)
# save the initial z as reference
S.lab.z0=z0
# assign constant velocity to the wall
S.lab.topWall.vel=(0,0,-.05)
# add to the simulation
# when the wall goes down to z=1, make it move upwards again
elif S.lab.stage==1 and S.lab.topWall.pos[2]<1:
S.lab.stage+=1
S.lab.topWall.vel*=-1
# move until the initial position is reached again
elif S.lab.stage==2 and S.lab.topWall.f[2]==0:
S.lab.stage+=1
# stop Wall velocity
S.lab.topWall.vel=(0,0,0)
# print z0 after compression
print('z0' , S.lab.topWall.pos[2])
# remove Wall
S.dem.par.remove(S.lab.topWall.id)
# command to start the ballast incertion
elif S.lab.stage==3 and S.lab.inletSL.dead and woo.utils.unbalancedForce(S)<0.05:
S.lab.stage+=1
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
z1=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
# place Wall just above the top
S.lab.topWall1=Wall.make(z1,axis=2,sense=-1,mat=S.lab.wallMat)
# print z1 before compression
print('z1' , z1)
# save the initial z as reference
S.lab.z1=z1
# assign constant velocity to the wall
S.lab.topWall1.vel=(0,0,-.05)
# add to the simulation
# when the wall goes down to z=1, make it move upwards again
elif S.lab.stage==4 and S.lab.topWall1.pos[2]<2:
S.lab.stage+=1
S.lab.topWall1.vel*=-1
# move until the initial position is reached again
elif S.lab.stage==5 and S.lab.topWall1.f[2]==0:
# everything done, nothing will happen with stage==3
S.lab.stage+=1
# stop Wall velocity
S.lab.topWall1.vel=(0,0,0)
# print z1 after compression
print('z1' , S.lab.topWall1.pos[2])
# remove Wall
S.dem.par.remove(S.lab.topWall1.id)
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
z2=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
print('z2' , z2)
# create the box
zBox=z2+zdor*.7
vx = (xmaxbox/2.)-(xdor/2.)
vy = (ymaxbox/2.)-(ydor/2.)
# for n in nn:
# lock Y to fix the dormente (.62 , .22 , .15)
# n.pos[2]==zBox make oscillation starts online after box reaches z2
# if abs(n.pos[2])<=zBox:
# n.dem.impose=HarmonicOscillation(freq=100.,amp=.02,dir=(0,0,1))
# print(n.pos[2])
for p in pp:
if isinstance(p.shape,Facet): p.shape.visible=False
else: p.shape.wire=False
S.lab.collider.noBoundOk=True
# S.lab.collider.sortThenCollide=True
for n in nn: DemData.setOriMassInertia(n)
S.lab.ncycle=S.step
print('ncycle' , S.lab.ncycle)
elif S.lab.stage==6 and S.step>=(S.lab.ncycle+2000):
S.lab.stage+=1
z3=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
print('z3' , z3)
# point of Box insertion
zBox3=z3+zdor+0.05
# create the box
# add to the simulation
S.lab.boxNode=Node(pos=(1,1,zBox3),dem=ClumpData(blocked='xyzXYZ'))
# assign velocity to the box
S.lab.boxNode.dem.impose=HarmonicOscillation(freq=100.,amp=.02,dir=(0,0,1))



commented Jan 17 by (46,390 points)

Hi Alan, could you send me the dormente.nf file as well? You can attach here, if you don't mind it being publicly accessible, or by e-mail? It will be better for me to see what you're after :) Thanks!

commented Jan 17 by (200 points)

Hi. I couldn't paste here, what is your e-mail? Thanks

commented Jan 18 by (46,390 points)

vaclav.smilauer@woodem.eu

commented Jan 18 by (200 points)

Send you the e-mail, just confirm to me if it arrived. Thanks

answered Jan 24 by (46,390 points)
Hi, I received the e-mail, thank you. There is indeed the explosion, so a few points to mention.

1. The explosion occurs because timestep changes abruptly by about 3 orders of magnitude. Since DynDt runs every 100 steps by default, one step is enough to lose it. The best remedy is, just after adding the mesh, say something like this: S.dynDt.stepPeriod=1 (run every step), S.dt*=1e-4 (reduce timestep right away). It will not (seems like, explosion, but "slow") work anyway, there might be some problem with your mesh which I did not check really.

2. Why do you put the dormente on the top and yet another bounding box on the top of it? The dormente has about the size of 2-3 spheres, so you could just as well impose force on those few spheres. In another words, the dormente is disproportionally small, thus the timestep drops heavily and your simulation will take forever.

HTH, Vaclav
commented Jan 30 by (200 points)
1. I adjusted the dtSafety to 0.05, and fixed some problems at the nmesh, so it could make contact with the spheres and removed the box.

2. Well, I was trying to represent something like this:

But now I using only the dormant (nmesh) and applying a Harmonic Oscillation on it. Also, this is a generic model, and I was just trying to make contact between the dormand and the spheres, but later I'll try with some smallers spheres.

The code is now at this stage, I could check that there is contact between the mesh and the spheres (at Display). But I coudn't make the mesh fall and stop at the top of the spheres, so I add it right above the spheres and apply the load.

Anyway, thank you for your help!

from woo.core import *
from woo.dem import *
from minieigen import *
from woo.fem import *
import woo.triangulated
from math import pi
import woo.models, woo

woo.master.usesApi=10102

# ioneMask=2 >> Particle groups which have bits in loneMask in common (i.e. (A.mask & B.mask & loneMask)!=0) will not have contacts between themselves

# contact model plus material for particles
model=woo.models.ContactModelSelector(name='linear',damping=.4,mats=[FrictMat(young=5e5,ktDivKn=.2,tanPhi=.4)])
# same material for walls as for particles, but without friction
S.lab.wallMat=model.mats[0].deepcopy(tanPhi=0.)

# intial box size
xmaxbox = 1.
ymaxbox = .5
zmaxbox = 1.
box=AlignedBox3((0,0,0),(xmaxbox,ymaxbox,zmaxbox))

# subballast size
slx0 = 0.
sly0 = 0.
slz0 = 0.
slxmax = 1.
slymax = .5
slzmax = .5
slastro=AlignedBox3((slx0,sly0,slz0),(slxmax,slymax,slzmax))

# ballast size
lx0 = 0.
ly0 = 0.
lz0 = slzmax*.75
lxmax = slxmax
lymax = slymax
lzmax = slzmax*1.5
lastro=AlignedBox3((lx0,ly0,lz0),(lxmax,lymax,lzmax))

# use something else than spheres here, if you like
gen2=PsdSphereGenerator(psdPts=[(.05,0),(.075,.33),(.1,.67),(.125,1)])
gen=PsdSphereGenerator(psdPts=[(.05,0),(.075,.5),(.1,1)])

# every 100 steps
ctrlStep=100

# dormant size
xdor = .62
ydor = .22
zdor = .15

# dormant material
m1 = woo.utils.defaultMaterial()
m1.young = 2.6E9
m1.density = 700 # make it faster

# wall box without the top (will be added later)

S.engines=DemField.minimalEngines(model=model)+[
# generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
# check what's happening periodically
PyRunner(ctrlStep,'checkProgress(S)')]+[IntraForce([In2_Tet4_ElastMat(),In2_Facet(),In2_Membrane_FrictMat(bending=True)])]

S.lab.stage=0
# avoid warning when we change the value
S.lab._setWritable('stage')
# things that we want to look at during the simualtion
S.plot.plots={'dz':('Fz0','Fz1')}

S.saveTmp()

# callback function from PyRunner
def checkProgress(S):
# if the top wall is there, save its position and acting force
if S.lab.stage in (1,2):
# graphic adjustment z0 x Fz0
elif S.lab.stage in (4,5):
# graphic adjustment z1 x Fz1
if S.lab.stage==0 and S.lab.inletL.dead and woo.utils.unbalancedForce(S)<0.05:
S.lab.stage+=1
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
z0=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
# place Wall just above the top
S.lab.topWall=Wall.make(z0,axis=2,sense=-1,mat=S.lab.wallMat)
# print z0 before compression
print('z0' , z0)
# save the initial z as reference
S.lab.z0=z0
# assign constant velocity to the wall
S.lab.topWall.vel=(0,0,-.05)
# add to the simulation
# when the wall goes down to z=1, make it move upwards again
elif S.lab.stage==1 and S.lab.topWall.pos[2]<(slzmax*.5):
S.lab.stage+=1
S.lab.topWall.vel*=-1
# move until the initial position is reached again
elif S.lab.stage==2 and S.lab.topWall.f[2]==0:
S.lab.stage+=1
# stop Wall velocity
S.lab.topWall.vel=(0,0,0)
# print z0 after compression
print('z0' , S.lab.topWall.pos[2])
# remove Wall
S.dem.par.remove(S.lab.topWall.id)
# command to start the ballast incertion
elif S.lab.stage==3 and S.lab.inletSL.dead and woo.utils.unbalancedForce(S)<0.05:
S.lab.stage+=1
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
z1=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
# place Wall just above the top
S.lab.topWall1=Wall.make(z1,axis=2,sense=-1,mat=S.lab.wallMat)
# print z1 before compression
print('z1' , z1)
# save the initial z as reference
S.lab.z1=z1
# assign constant velocity to the wall
S.lab.topWall1.vel=(0,0,-.05)
# add to the simulation
# when the wall goes down to z=1, make it move upwards again
elif S.lab.stage==4 and S.lab.topWall1.pos[2]<(lzmax*.5):
S.lab.stage+=1
S.lab.topWall1.vel*=-1
# move until the initial position is reached again
elif S.lab.stage==5 and S.lab.topWall1.f[2]==0:
# everything done, nothing will happen with stage==3
S.lab.stage+=1
# stop Wall velocity
S.lab.topWall1.vel=(0,0,0)
# print z1 after compression
print('z1' , S.lab.topWall1.pos[2])
# remove Wall
S.dem.par.remove(S.lab.topWall1.id)
# find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
S.lab.z2=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
print('z2' , S.lab.z2)
# create the box
zBox=S.lab.z2+zdor*.05
vx = (xmaxbox/2.)-(xdor/2.)
vy = (ymaxbox/2.)-(ydor/2.)
# S.dynDt.stepPeriod=1
# S.dt*=1e-4
S.lab.ncycle=S.step+5000
print('start at ncycle' , S.lab.ncycle)
for p in pp:
if isinstance(p.shape,Facet): p.shape.visible=False
else: p.shape.wire=False
S.lab.collider.noBoundOk=True
for n in nn: DemData.setOriMassInertia(n)
# S.lab.ncycle=S.step
# print('ncycle' , S.lab.ncycle)
S.lab.nn = nn
elif S.lab.stage==6 and S.step>=(S.lab.ncycle):
S.lab.stage+=1

for n in S.lab.nn:
if abs(n.pos[2])>=S.lab.z2+(.8*zdor):
# print('start at n.pos' , n.pos[2])
n.dem.impose=HarmonicOscillation(freq=100.,amp=.02,dir=(0,0,-1))