yade-users team mailing list archive

[Alexandre Sac-Morane <question707162@xxxxxxxxxxxxxxxxxxxxx>]

New question #707162 on Yade:
https://answers.launchpad.net/yade/+question/707162

Dear all,

I am trying to use PotentialBlocks to model particles with nonregular shapes. In the following script, I define 2 PotentialBlocks from a sphere. I know there is a yade function, but it is a first example. The idea is then to have a nonregular shape. One grain is fixed, and the other is falling to it. The grains are defined with an angular discretization ('n_theta') related to the number of planes considered.

It appears when this discretization is equal to 7, the code is working well. But when this discretization is equal to at least 8, a Segmentation fault occurs. Grains are correctly generated but when the contact occurs the simulation stops.

What do you think about it? Is it because of the computational cost (increasing the number of planes)?

Thanks in advance for your answer.
Regards
Alexandre Sac-Morane

Here, my mwe :


from yade import utils
from potential_utils import *
import math

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Material Parameters
m = O.materials.append(FrictMat(young=-1, poisson=-1, frictionAngle=atan(0.5), density=2000, label='frictMat'))

Kn = 1e7  #Pa/m
Ks = Kn * 2 / 3  #Pa/m
kn_i = 5 * Kn
ks_i = 5 * Ks

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([PotentialBlock2AABB()], verletDist=0.00),
        InteractionLoop(
                [Ig2_PB_PB_ScGeom(calContactArea=True)],
                [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=kn_i, ks_i=ks_i, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
                [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', initialOverlapDistance = 1e-6, allowViscousAttraction=False)]
        ),
        NewtonIntegrator(damping=0.0, exactAsphericalRot=True, gravity=[0, 0, -9.81], label='newton')
]

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Parameters
radius = 1 # the radius of grains
r = 0.1
n_theta = 7 #7 max with this definition on my laptop
n_phi = int(n_theta/2+1)

# creation of a discretized sphere
a_L = []
b_L = []
c_L = []
d_L = []
for i_phi in range(n_phi):
	phi = math.pi*i_phi/(n_phi-1)
	# top plane
	if i_phi == 0:
		a_L.append(0) # x coordinate of the norm of the plane
		b_L.append(0) # y coordinate of the norm of the plane
		c_L.append(1) # z coordinate of the norm of the plane
		d_L.append(radius-r) # distance to the center of the plane
	# bottom plane
	elif i_phi == n_phi-1:
		a_L.append( 0) # x coordinate of the norm of the plane
		b_L.append( 0) # y coordinate of the norm of the plane
		c_L.append(-1) # z coordinate of the norm of the plane
		d_L.append(radius-r) # distance to the center of the plane
	# current plane
	else :
		if int(n_theta*math.sin(phi)) != 0 :
			for i_theta in range(int(n_theta*math.sin(phi))):
				theta = 2*math.pi*i_theta/int(n_theta*math.sin(phi))
				a_L.append(math.cos(theta)*math.sin(phi)) # x coordinate of the norm of the plane
				b_L.append(math.sin(theta)*math.sin(phi)) # y coordinate of the norm of the plane
				c_L.append(                math.cos(phi)) # z coordinate of the norm of the plane
				d_L.append(radius-r) # distance to the center of the plane

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# potential_utils.potentialblock: Generic Potential Block defined by a,b,c,d

# grain 1
g1 = Body()
g1.aspherical = True
g1.shape = PotentialBlock(
        a=a_L,
        b=b_L,
        c=c_L,
        d=d_L,
        r=r,
        R=0.0,
        AabbMinMax=True,
        isBoundary=True
)
utils._commonBodySetup(g1, g1.shape.volume, g1.shape.inertia, material='frictMat', pos=[0, 0, 0], fixed=True)
g1.state.pos = [0, 0 , -radius]
g1.state.ori = g1.shape.orientation
O.bodies.append(g1)

# grain 2
g2 = Body()
g2.aspherical = True
g2.shape = PotentialBlock(
        a=a_L,
        b=b_L,
        c=c_L,
        d=d_L,
        r=r,
        R=0.0,
        AabbMinMax=True,
        isBoundary=False
        )
utils._commonBodySetup(g2, g2.shape.volume, g2.shape.inertia, material='frictMat', pos=[0, 0, 0], fixed=False)
g2.state.pos = [0, 0 , 1.1*radius]
g2.state.ori = g2.shape.orientation
O.bodies.append(g2)

# ------------------------------------------------------------------------------------------------------------------------------------------ #
# Timestep
O.dt = 5e-5

O.run()


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.