← Back to team overview

yade-users team mailing list archive

Re: [Question #695558]: Dynamics of pfacet object does not appear to be correct

 

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

Description changed to:
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

# ---------------------------------------------------------------------------------------------------------------- code
from yade import plot, utils
from yade.gridpfacet import *

# ---------------------------------------------------------------------- defining material
bristle_radius = 1e-3
bristle_length = 0.1
bristle_seg_no = 5

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

grid_radius = 1e-3
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(),
    ]),
    InteractionLoop(
        [
            Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
            Ig2_GridNode_GridNode_GridNodeGeom6D(),
            Ig2_PFacet_PFacet_ScGeom(),
            Ig2_GridConnection_PFacet_ScGeom(),
            Ig2_Sphere_PFacet_ScGridCoGeom(),
        ],
        [
            Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False),
            Ip2_FrictMat_FrictMat_FrictPhys()
        ],
        [
            Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
            Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
        ]
    ),
    NewtonIntegrator(gravity = (0,0,-10),damping = 0.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 + origin[2]],
    [-x_len,  0,     0 + origin[2]],
    [0.0, y_len, z_len + origin[2]],
    [0.0,-y_len, z_len + origin[2]],
]

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


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
# ---------------------------------------------------------------------- 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.