Two material

Hi,

I working on compaction, with two different material. The material type is the same, the difference is the radius and the density.
My simulation works with single materials, however when use two material the particle began to fly away during gravity filling.

Best regards,

MWE:

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from yade import pack, export
import os
readParamsFromTable(comp_press=0.81e8,h_tab=2.02,m_tab=0.2117, r_tab=5.015,wCCS=0.02,tab_porosity=15,tab_height=1,save=1)
from yade.params.table import *

o = Omega()
save=save
# Physical parameters
fr = 0.41
rho_PH101 = 1561
rho_CCS =1403
D_PH101 = 7.9e-5
r1_PH101 = D_PH101/2
D_CCS = 5.4e-5
r1_CCS = D_CCS/2
#r2 = Diameter/2
k1 = 10000
kp = 140000
kc = k1 * 0.1
ks = k1 * 0.1
Chi1 = 0.34

o.dt = 1.0e-8

PhiF1=0.999

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

Comp_press_up= comp_press
Comp_force_up=Comp_press_up*cross_area

Comp_press_lp= comp_press
Comp_force_lp=Comp_press_lp*cross_area

# Add material
matPH101 = O.materials.append(LudingMat(frictionAngle=fr, density=rho_PH101, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))
matCCS = O.materials.append(LudingMat(frictionAngle=fr, density=rho_CCS, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp_PH101=pack.SpherePack()
sp_PH101.makeCloud((-8*D_PH101,-8*D_PH101,-50*D_PH101),(8*D_PH101,8*D_PH101,40.0*D_PH101), rMean=r1_PH101,rRelFuzz=0.18,num=round(10134))
sp_PH101.toSimulation(material=matPH101,color=(0,1,1))
sp_CCS=pack.SpherePack()
sp_CCS.makeCloud((-8*D_PH101,-8*D_PH101,-50*D_PH101),(8*D_PH101,8*D_PH101,40.0*D_PH101), rMean=r1_CCS,rRelFuzz=0.15,num=round(740))
sp_CCS.toSimulation(material=matCCS,color=(1,0,1))

walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=matPH101))

# 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,-9.81]),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]

def checkForce():
    if O.iter < 3000000:
        return
    if unbalancedForce() > 1:
        return
    global upper_punch
    upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for i in upper_punch:
        body= O.bodies[i]
        body.state.vel = (0,0,-0.02)
    global lower_punch
    lower_punch= O.bodies.append(geom.facetCylinder((0,0,(-Cyl_height/2)-0.0001),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for n in lower_punch:
        body= O.bodies[n]
        body.state.vel = (0,0,0.02)
    O.engines = O.engines + [PyRunner(command='storeData()', iterPeriod=25000)]+ [PyRunner(command='saveData()', iterPeriod=100000)]
    fCheck.command = 'unloadPlate()'