yade-users team mailing list archive

[Rohit John <question695558@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello,

I am simulating the interaction between a brush (bristles made of grid nodes and grid connections) and an object (made of pfacet). The bristles are arranged to form a square at its  base. The objects takes the form of a wedge placed symmetrically between the bristles. The object is made to move into the brush

I expect the object to bounce back up, but it is starting to rotate when it contacts the brush. Since the system is symmetric, I do not expect this rotation to happen. I could not find the bug in my code. I suppose there could be something in the dynamics I am missing.

Kind regards,
Rohit John
The minimum working code
# ---------------------------------------------------------------------------------------------------------------------------------------
from yade import plot, utils
from yade.gridpfacet import *

# user define function
# ---------------------------------------------------------------------------------------------- add motion engines
def add_motion_engines(ids, linear_velocity, angular_velocity, center_of_rotation):
    '''
    Adds the engines which sets the initial velocities of a target. This is done because for some reason we cannot graph the velocities for non spherical objects by manually setting the velocities
    '''
    linear_velocity_magnitude = L2_norm(linear_velocity)
    linear_velocity_vector    = normalise(linear_velocity)

    angular_velocity_magnitude = L2_norm(angular_velocity)
    angular_velocity_vector    = normalise(angular_velocity)

    O.engines += [
        CombinedKinematicEngine(
            ids   = ids,
            label ='combined_motion_engine') + 
            
        TranslationEngine(
            translationAxis = linear_velocity_vector,
            velocity        = linear_velocity_magnitude,
            label           = 'translation_engine') + 
                
        RotationEngine(
            rotationAxis     = angular_velocity_vector, 
            angularVelocity  = angular_velocity_magnitude, 
            rotateAroundZero = True, 
            zeroPoint        = center_of_rotation,
            label            = 'rotation_engine'
            ),

        PyRunner(command = "motion_engine_stopper()", iterPeriod = 1, label = "stopper")
        
    ]

# ---------------------------------------------------------------------------------------------- motion_engine_stopper
def motion_engine_stopper():
    combined_motion_engine.dead = True
    translation_engine.dead = True
    rotation_engine.dead    = True
    stopper.iterPeriod      = int(1e14)
    print("Engines stopped")

# ---------------------------------------------------------------------------------------------- L2_norm
def L2_norm(vector):
    mag = 0
    for i in vector:
        mag += i**2

    return mag**0.5

# ---------------------------------------------------------------------------------------------- normalise
def normalise(vector):
    mag = L2_norm(vector)
    if mag == 0:
        return [0,0,0]
    res = [i/mag for i in vector]
    return res

# ---------------------------------------------------------------------- defining material
bristle_radius = 1e-3
bristle_length = 0.1
brush_x_dim    = 0.1
brush_y_dim    = 0.02
brush_position = [0,0.0025,0]
brush_ori      = Quaternion((1,0,0), 0)

bristle_x_no   = 10
bristle_y_no   = 2
bristle_seg_no = 5
bristle_x_density = bristle_x_no / brush_x_dim
bristle_y_density = bristle_y_no / brush_y_dim

bristle_young    = 3e9
bristle_poisson  = 0.3
bristle_density  = 1000
bristle_friction = radians(30)

covar = 0
bristle_tip_spread_covariance = [[covar, 0], [0, covar]]

target_dimension = [0.1,0.1,0.1]
target_position  = [0,0,0]
grid_radius     = 1e-3
target_ori       = Quaternion((1,0,0), pi/4)
target_young    = 30e9
target_poisson  = 0.3
target_density  = 1000
target_friction = radians(30)

# ---------------------------------------------------------------------- material
O.materials.append(
    FrictMat(
        young   = bristle_young,
        poisson = bristle_poisson,
        density = bristle_density,
        label   = 'cyl_ext_mat',
        frictionAngle = bristle_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young   = bristle_young,
        poisson = bristle_poisson,
        density = bristle_density,
        label   = 'cyl_int_mat',

        frictionAngle     = bristle_friction,
        momentRotationLaw = True,
        normalCohesion    = 1e40,
        shearCohesion     = 1e40,
    )
)

O.materials.append(
    FrictMat(
        young   = target_young,
        poisson = target_poisson,
        density = target_density,
        label   = 'pfacet_ext_mat',
        frictionAngle = target_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young   = 100*target_young,
        poisson = target_poisson,
        density = target_density,
        label   = 'pfacet_int_mat',

        frictionAngle     = target_friction,
        momentRotationLaw = True,
        normalCohesion    = 1e40,
        shearCohesion     = 1e40,
    )
)

# ---------------------------------------------------------------------- Engines
O.engines = [
        ForceResetter(),

        InsertionSortCollider([
            Bo1_GridConnection_Aabb(),
            Bo1_PFacet_Aabb(),
            Bo1_Sphere_Aabb(),
        ]),

        InteractionLoop(
            [
                Ig2_PFacet_PFacet_ScGeom(),
                Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
                Ig2_GridNode_GridNode_GridNodeGeom6D(),
                Ig2_GridConnection_PFacet_ScGeom(),
                Ig2_Sphere_PFacet_ScGridCoGeom(),
                Ig2_Sphere_Sphere_ScGeom()
            ],
            [
                Ip2_FrictMat_FrictMat_FrictPhys(),
                Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(
                    setCohesionNow = True, 
                    setCohesionOnNewContacts = False
                    ),
            ],
            [
                Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
                Law2_ScGeom_FrictPhys_CundallStrack(),
                Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
                Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
            ],
        ),
        NewtonIntegrator(gravity = [0.0,0.0,0.0], damping = 0)
    ]   
# ---------------------------------------------------------------------- bodies 
# custome brush
nodesIds=[]
cylIds=[]

# bristle root location
base_pos = [
    [0.003, 0.003, 0],
    [0.003, -0.003, 0],
    [-0.003, 0.003, 0],
    [-0.003, -0.003, 0],
    ]
dz = bristle_length / bristle_seg_no # grid connection length

for i in base_pos:
    bristle = []
    for j in range(bristle_seg_no + 1):
        bristle.append([
            i[0],
            i[1],
            bristle_length - (i[2] + dz * j) 
        ])
    cylinderConnection(
        bristle,
        bristle_radius,
        nodesIds,
        cylIds,
        color=[1,0,0],
        intMaterial='cyl_int_mat',
        extMaterial='cyl_ext_mat'
    )
    O.bodies[nodesIds[-1]].state.blockedDOFs='xyzXYZ'
        
# manaul target
origin = [
    0,
    0,
    bristle_length + bristle_radius, + grid_radius
]

x_len = 0.03
y_len = 0.05
z_len = 0.05
nodes = [
    [x_len,  0, 0 + bristle_length + bristle_radius + grid_radius],
    [-x_len, 0, 0 + bristle_length + bristle_radius + grid_radius],
    [0.0, y_len, z_len + bristle_length + bristle_radius + grid_radius],
    [0.0,-y_len, z_len + bristle_length + bristle_radius + grid_radius]
]

target_pos = [0,0,0] # centroid

for j in range(3):
    for i in nodes:
        target_pos[j] += i[j]
    target_pos[j] / 4
    

pfacet_nodes = [
    [nodes[0], nodes[1], nodes[2]],
    [nodes[0], nodes[1], nodes[3]]
]

pnode  = []
pcyl   = []
pfacet = []
for i in pfacet_nodes:
    print(i)
    pfacetCreator1(
        i,
        grid_radius, 
        nodesIds = pnode, 
        cylIds   = pcyl, 
        pfIds    = pfacet, 
        wire     = False, 
        fixed    = False, 
        color    = [0.5,0.5,0.5],
        materialNodes = 'pfacet_int_mat', 
        material      = 'pfacet_ext_mat', 
        )

target_ids = pnode + pcyl + pfacet
# # ---------------------------------------------------------------------- motion engine
target_lin_vel = [0,0,-1]
target_ang_vel = [0,0,0]

add_motion_engines(
    ids = target_ids,
    linear_velocity  = target_lin_vel,
    angular_velocity = target_ang_vel,
    center_of_rotation = target_pos
)

# ----------------------------------------------------------------------------- graphing and additional engines
graph_iter = 1000
plot.plots = {
    't':(
        'wx_avg', 'wy_avg', 'wz_avg',
        )}

# ----------------- initialising the variables
# These variables will be added each iteration and its average will be taken (dividing my the number of iteration) when the graph is plotted
wx_avg = 0
wy_avg = 0
wz_avg = 0
# bodyList = target_body.node_id_list
bodyList = target_ids

def get_avg_state():    
    global wx_avg, wy_avg, wz_avg

    wx_temp = 0 
    wy_temp = 0 
    wz_temp = 0 

    for i in bodyList:
        wx_temp += O.bodies[i].state.angVel[0]
        wy_temp += O.bodies[i].state.angVel[1]
        wz_temp += O.bodies[i].state.angVel[2]
        
    # averaging for all bodies
    no_of_bodies = len(bodyList)

    wx_temp /= no_of_bodies
    wy_temp /= no_of_bodies
    wz_temp /= no_of_bodies    

    # addind to the cumulative states. This will later be divided by graph_iter to get the average

    wx_avg += wx_temp
    wy_avg += wy_temp
    wz_avg += wz_temp

def graph():
    global wx_avg, wy_avg, wz_avg

    wx_avg /= graph_iter
    wy_avg /= graph_iter
    wz_avg /= graph_iter

    plot.addData(t = O.time,  
        wx_avg = wx_avg, wy_avg = wy_avg, wz_avg = wz_avg,
        )

plot.plot()
# ---------------------------------------------------------------------- plotting and simulation stop engines
end_time = 4e-1
O.engines += [
    PyRunner(command = 'graph()', iterPeriod = graph_iter),
    PyRunner(command = 'get_avg_state()', iterPeriod = 1),
]

O.dt = 5e-7
O.saveTmp()

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