yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #25369
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
Status: Answered => Open
Rohit John is still having a problem:
Dear Bruno,
I have found one more bug. The angular momentum value of the clump
returns (0,0,0) even though inertia is not 0 and angular velocity is not
0. This happens when I set the mass of the gridConnections to 0.
However, if I do not do this, I get a non zero value for the angular
momentum.
In the simulation, two sphere moving the opposite direction. They are
not collinear, they are offset. One is offset in the +y while the other
in the -y direction. The are set to contact a pfacet cube at the same
moment of time on opposite faces. I wanted to check if the angular
momentum was conserved. It turns out that it is not conserved.
I am using ubuntu 18.04. I found this error in yade-2021-03-22.git-096f9a3, yade-2020-11-05.git-7d0ae90
and Yade 2018.02b. All of them showed this bug
Kind regards,
Rohit K. John
#---------------------------------------------------- GTS FILE. File name = 'cube.gts'
14 36 24 GtsSurface GtsFace GtsEdge GtsVertex
0.5 0.5 0.5
0.5 0.5 -0.5
0.5 -0.5 0.5
0.5 -0.5 -0.5
-0.5 0.5 0.5
-0.5 0.5 -0.5
-0.5 -0.5 0.5
-0.5 -0.5 -0.5
0.5 0.0 0.0
0.0 -0.5 0.0
0.0 0.0 0.5
-0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 -0.5
6 8
2 6
1 2
8 7
3 4
5 6
3 7
1 3
8 4
7 5
5 1
4 2
9 3
2 9
4 9
1 9
10 7
4 10
8 10
3 10
11 3
5 11
7 11
1 11
12 6
7 12
8 12
5 12
13 2
5 13
6 13
1 13
14 8
2 14
6 14
4 14
7 21 23
7 17 20
6 25 28
9 33 36
8 13 16
3 29 32
12 14 15
5 15 13
3 16 14
9 18 19
4 19 17
5 20 18
11 22 24
10 23 22
8 24 21
4 26 27
1 27 25
10 28 26
6 30 31
2 31 29
11 32 30
2 34 35
1 35 33
12 36 34
#------------------------------------------------------------------------------------------- YADE Script
from yade.gridpfacet import *
from yade import geom, utils
from yade import plot
import sys, os
sys.path.append(".")
# ---------------------------------------------------------------------------- input parameter
# ----------------------------------------------------- target
target_young = 50e9
target_density = 1000
target_poisson = 0.3
target_friction = radians(30)
p_radius = 5e-2
# ---------------------------------------------------------------------------------------------------------- Materials
target_int_mat = 'pfacet_int_mat'
target_ext_mat = 'pfacet_ext_mat'
O.materials.append(
FrictMat(
young = target_young,
poisson = target_poisson,
density = target_density,
label = target_ext_mat,
frictionAngle = target_friction,
)
)
O.materials.append(
CohFrictMat(
young = target_young,
poisson = target_poisson,
density = target_density,
label = target_int_mat,
frictionAngle = target_friction,
normalCohesion = 3e100,
shearCohesion = 3e100,
momentRotationLaw = True,
)
)
# ---------------------------------------------------------------------------------------------------------- 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(),
],
[
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(),
],
)
]
# ---------------------------------------------------------------------------------------------------------- objects
# ---------------------------------------------------------------------- target
(
pnode,
pcyl,
pfacet
) = gtsPFacet(
'cube.gts',
radius = p_radius,
shift = (0,0,0),
scale = 1,
wire = False,
fixed = False,
color = [0.1,0.5,0.1],
materialNodes = 'pfacet_int_mat',
material = 'pfacet_ext_mat',
)
target_ids = pnode + pcyl + pfacet
for i in pcyl:
O.bodies[i].state.mass = 0
target_clump_ID = O.bodies.clump(target_ids)
O.bodies[target_clump_ID].state.inertia = Vector3([0.01,0.01,0.01])
O.bodies[target_clump_ID].state.mass = 1
print(target_clump_ID)
# ---------------------------------------------------------------------- spheres
sp1 = sphere([ 0.7, 0.25, 0], 5e-2, material = target_ext_mat)
sp2 = sphere([-0.7, -0.25, 0], 5e-2, material = target_ext_mat)
sp1_ID = O.bodies.append(sp1)
sp2_ID = O.bodies.append(sp2)
O.bodies[sp1_ID].state.vel = [-1,0,0]
O.bodies[sp2_ID].state.vel = [ 1,0,0]
# ----------------------------------------------------------------------------- Additional engines
ids = target_clump_ID
O.engines += [
NewtonIntegrator(gravity = [0,0,0], damping = 0),
PyRunner(command = 'graph()', iterPeriod = 1000)
]
# ---------------------------------------------------------------------- plotting
plot.plots = {'t':('wx', 'wy', 'wz'), 't1':('Lx', 'Ly', 'Lz')}
def graph():
w = O.bodies[ids].state.angVel
pos_sp1 = O.bodies[sp1_ID].state.pos
vel_sp1 = O.bodies[sp1_ID].state.vel
angMom_sp1 = O.bodies[sp1_ID].state.angMom
mass_sp1 = O.bodies[sp1_ID].state.mass
L_sp1 = angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1)
pos_sp2 = O.bodies[sp2_ID].state.pos
vel_sp2 = O.bodies[sp2_ID].state.vel
angMom_sp2 = O.bodies[sp2_ID].state.angMom
mass_sp2 = O.bodies[sp2_ID].state.mass
L_sp2 = angMom_sp2 + mass_sp2 * pos_sp2.cross(vel_sp2)
angVel_clump = O.bodies[target_clump_ID].state.angVel
mmoi_clump = O.bodies[target_clump_ID].state.inertia
mmoi_mat_clump = Matrix3([
mmoi_clump[0], 0, 0,
0, mmoi_clump[1], 0,
0, 0, mmoi_clump[2]
])
L_clump = mmoi_mat_clump * angVel_clump
L_tot = L_clump + L_sp1 + L_sp2
angMom_clump = O.bodies[target_clump_ID].state.angMom
print(L_tot)
plot.addData(
t = O.time, t1 = O.time,
wx = w[0], wy = w[1], wz = w[2],
Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2],
)
plot.plot()
# ----------------------------------------------------------------------------- Sim controller
O.dt = 2e-7#utils.PWaveTimeStep()
O.saveTmp()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
