Welcome to ask.woodem.org. You may post when you login through your GitHub account.

Contact Particle + Nmesh and Nmesh + triangulated.box

0 votes
17 views
asked Jan 17 by alanbds (200 points)

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

S=woo.master.scene=Scene(fields=[DemField(gravity=(0,0,-10),loneMask=0)],dtSafety=.5)
# 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
mask=0b001
mask2=0b011

# wall box without the top (will be added later)
S.dem.par.add(Wall.makeBox(box=box,which=(1,1,1,1,1,0),mat=S.lab.wallMat))

S.engines=DemField.minimalEngines(model=model)+[
	# generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
	BoxInlet(stepPeriod=ctrlStep,box=slastro,generator=gen2,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inletL'),
	BoxInlet(stepPeriod=ctrlStep,box=lastro,generator=gen,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inletSL',dead=True),
	# check what's happening periodically
	PyRunner(ctrlStep,'checkProgress(S)')]+[IntraForce([In2_Tet4_ElastMat(),In2_Facet(),In2_Membrane_FrictMat(bending=True)])]

# stage: 0=gravity deposition subballast, 1=oedometric loading, 2=oedometric unloading, 3=closing subballast, 4=gravity deposition ballast, 5=oedometric loading, 6=oedometric unloading, 7=box loading,
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
		S.plot.addData(dz=S.lab.z0-S.lab.topWall.pos[2],Fz0=S.lab.topWall.f[2],i=S.step,t=S.time)
	elif S.lab.stage in (4,5):
		# graphic adjustment z1 x Fz1
		S.plot.addData(dz=S.lab.z1-S.lab.topWall1.pos[2],Fz1=S.lab.topWall1.f[2],i=S.step,t=S.time)
	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
		S.dem.par.add(S.lab.topWall)
		# 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
		# reverse the loading sense
		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
		S.lab.inletSL.dead=False
	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
		S.dem.par.add(S.lab.topWall1)
		# 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
		# reverse the loading sense
		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.)
		nn,pp=woo.utils.importNmesh('dormente.nf',mat=m1,mask=mask,trsf=lambda v: v+Vector3(vx,vy,zBox),dem=S.dem,surfHalfThick=0.005)
		# nn,pp=woo.utils.importNmesh('dormente.nf',mat=m1,mask=mask,dem=S.dem)
		# 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
		S.lab.Box=woo.triangulated.box(dim=(.1,.1,.1),center=(1,1,zBox3),which=(1,1,1,1,1,1),fixed=True,mat=S.lab.wallMat,mask=mask,wire=False)
		# add to the simulation
		S.lab.boxNode=Node(pos=(1,1,zBox3),dem=ClumpData(blocked='xyzXYZ'))
		S.dem.par.addClumped(S.lab.Box,centralNode=S.lab.boxNode)
		# assign velocity to the box
		S.lab.boxNode.dem.impose=HarmonicOscillation(freq=100.,amp=.02,dir=(0,0,1))	
		

 

commented Jan 17 by eudoxos (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 alanbds (200 points)

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

commented Jan 18 by eudoxos (46,390 points)

vaclav.smilauer@woodem.eu

commented Jan 18 by alanbds (200 points)

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

1 Answer

0 votes
answered Jan 24 by eudoxos (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 alanbds (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

S=woo.master.scene=Scene(fields=[DemField(gravity=(0,0,-10),loneMask=0)],dtSafety=.05)
# 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
mask = 0b001
mask2 = 0b011

# wall box without the top (will be added later)
S.dem.par.add(Wall.makeBox(box=box,which=(1,1,1,1,1,0),mat=S.lab.wallMat))

S.engines=DemField.minimalEngines(model=model)+[
    # generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
    BoxInlet(stepPeriod=ctrlStep,box=slastro,generator=gen2,materials=[model.mats[0]],massRate=100,maxMass=0,maxAttempts=10000,atMaxAttempts='dead',label='inletL'),
    BoxInlet(stepPeriod=ctrlStep,box=lastro,generator=gen,materials=[model.mats[0]],massRate=100,maxMass=0,maxAttempts=10000,atMaxAttempts='dead',label='inletSL',dead=True),
    # check what's happening periodically
    PyRunner(ctrlStep,'checkProgress(S)')]+[IntraForce([In2_Tet4_ElastMat(),In2_Facet(),In2_Membrane_FrictMat(bending=True)])]

# stage: 0=gravity deposition subballast, 1=oedometric loading, 2=oedometric unloading, 3=closing subballast, 4=gravity deposition ballast, 5=oedometric loading, 6=oedometric unloading, 7=box loading,
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
        S.plot.addData(dz=S.lab.z0-S.lab.topWall.pos[2],Fz0=S.lab.topWall.f[2],i=S.step,t=S.time)
    elif S.lab.stage in (4,5):
        # graphic adjustment z1 x Fz1
        S.plot.addData(dz=S.lab.z1-S.lab.topWall1.pos[2],Fz1=S.lab.topWall1.f[2],i=S.step,t=S.time)
    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
        S.dem.par.add(S.lab.topWall)
        # 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
        # reverse the loading sense
        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
        S.lab.inletSL.dead=False
    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
        S.dem.par.add(S.lab.topWall1)
        # 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
        # reverse the loading sense
        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.)
        nn,pp=woo.utils.importNmesh('dormente.nf',mat=m1,mask=mask,trsf=lambda v: v+Vector3(vx,vy,zBox),dem=S.dem,surfHalfThick=0.005)
        # 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))
...