Increase particle radius

Asked by Mithushan Soundaranathan

Hi,

I am increasing my particle radius at each 100 iteration. But particle radius is not updating when using body.state.radius to update it. The change is not visible in the GUI.

Best,
Mithu

Here is my code:

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from scipy.integrate import odeint

o = Omega()

# Physical parameters
fr = 0.54
rho = 1050
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 100000
kp = 12.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.008
Cyl_height=0.045
cross_area=math.pi*(Tab_rad**2)

Comp_press= 0.5e8
Comp_force=Comp_press*cross_area
i=0

# swelling_mode function
def model(r,t,Q_max,rho_t,rho_w,r_0,D):
    Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
    drdt =((D*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
    return drdt

P=[1.45,rho,1000,r1,396.36e-12]
time_save=[]
radius_save=[]
#*************************************

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-18*Diameter),(3.0*Diameter,3.0*Diameter,18.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=527)
sp.toSimulation(material=mat1)
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, 0]),
  PyRunner(command='ParticleSwelling()', iterPeriod=100, label="fCheck"),
  #DeformControl(label="DefControl")
]

def ParticleSwelling():
    time_current=O.iter*o.dt
    time_save.append(time_current)
    radius_save_t=[]
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            t=[0,time_current-time_save[0]]
            r= odeint(model,P[3],t,args=(P[0],P[1],P[2],P[3],P[4],))
            b.state.radius = r[1]
            radius_save_t.append(b.state.radius)
    radius_save.append(radius_save_t)
    if time_current>8:
        O.pause()
        radius_data=pd.DataFrame(radius_save)
        radius_data.to_csv(r'compression_data_PH101.csv')

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hi,
use b.shape.radius (not state)
Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#2

Hi,

When I changed to b.shape.radius, I got following error message:

ArgumentError: Python argument types in
    None.None(Sphere, numpy.ndarray)
did not match C++ signature:
    None(yade::Sphere {lvalue}, double)

Best,
Mithu

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#3

Hi,

And one more thing, if update my particle radius in the terminal.

I am using following code:

for b in O.bodies:
 if isinstance(b.shape, Sphere):
    b.shape.radius=0.1

This code delete my spheres and add an additional cylinders (walls) with radius 0.1.

Revision history for this message
Jan Stránský (honzik) said :
#4

The first error says that it expects a number, but you provide a numpy.ndarray

> r= odeint(model,P[3],t,args=(P[0],P[1],P[2],P[3],P[4],))
> b.state.radius = r[1]

maybe here, r[1] might be something different than you want.

Concerning radii of cylinders, the condition "if isinstance(b.shape, Sphere)" is true also for cylinder, since Cylinder is inherited from Sphere [1].
Use
b.shape.__class__ == Sphere
instead to detect just spheres (and not derived classes too).

Cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#shape

Can you help with this problem?

Provide an answer of your own, or ask Mithushan Soundaranathan for more information if necessary.

To post a message you must log in.