Two contact law

Hello,

please provide a MWE [1] = full script, not just engines

> o.engines = [

please provide actual code.
o.engines should result in something like "'o' is not defined"

> do I need to add two different engine?

depends on what you mean.
There is exactly one O.engines for a simulation.
If you mean law functors [2], then yes, for two contact laws, you would need two Law functors

> I want to use Luding contact model between particles and The Hertz Mindlin non-slip contact model for particle-wall interactions. How can I use these different contact model in one simulation?

You want to use , Hertz Mindlin non-slip contact model, so you need:
- FrictMat
- Ip2_FrictMat_FrictMat_MindlinPhys
- Law2_ScGeom_MindlinPhys_Mindlin (?)

You want to use Luding contact model, so you need:
- LudingMat
- Ip2_LudingMat_LudingMat_LudingPhys
- Law2_ScGeom_LudingPhys_Basic

Normally it is OK to use different materials, Ips2 and Laws2 and their combinations.
E.g. I used to use Cpm contact model for sphere-sphere (spheres had Cpm material) and Cundal-Strack model for sphere-facets (facets had "only" FrictMat)

Unfortunately for you, the design of Yade and implementation decisions (the class hierarchy of materials mainly) made it impossible for one particle to have both Mindlin and Luding interactions.
Or it seems to me after looking briefly at the source code.
To achieve this, you would need some modifications of the source code...

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/user.html#functors-choice

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

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test

from yade import utils, plot, timing
from yade import pack
import pandas as pd

o = Omega()

# Physical parameters
fr = 0.52
rho = 1564 #kg/m3
Diameter = 0.011 #110 micrometer
r1 = Diameter
r2 = Diameter
k1 = 1005.0
kp = 10.0*k1
kc = k1 * 0.0
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-5

particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1

PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

#*************************************

# 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

sp=pack.SpherePack()
sp.makeCloud((-4.0*Diameter,-4.0*Diameter,-2.5*Diameter),(3.0*Diameter,3.0*Diameter,15.0*Diameter), rMean=Diameter/2.0, num=300)
sp.toSimulation()

##No of particles
count = 0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        count +=1
##Mass of particles
n_p=count

######################################################################
walls=O.bodies.append(yade.geom.facetCylinder((0,0,-0.02),radius=0.055,height=0.1,segmentsNumber=20,wallMask=6,material=mat1))

##Top punch##
#find max Z coordinate of your cloud
max_z = aabbExtrema()[1][2]
min_z = aabbExtrema()[0][2]

Top_punch=O.bodies.append(wall((0,0,max_z), axis=2, sense=-1, material=mat1))
O.bodies[Top_punch].state.vel = (0, 0, -0.01)

df = pd.DataFrame(columns=['time','Porosity','Maxoverlap','Stress'])
df.to_csv('PH102_Haustein_data.csv')
from csv import writer
def append_list_as_row(file_name, list_of_elem):
    # Open file in append mode
    with open(file_name, 'a+', newline='') as write_obj:
        # Create a writer object from csv module
        csv_writer = writer(write_obj)
        # Add contents of list as last row in the csv file
        csv_writer.writerow(list_of_elem)

# 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]),
  #VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=10000),
  #PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #PyRunner(command='print(voxelPorosity(resolution=600,start=aabbExtrema()[0],end=aabbExtrema()[1]))', realPeriod=15),
  PyRunner(command='checkdata()', realPeriod=1),
  DeformControl(label="DefControl")
]

def checkdata():
 if o.realtime<40:
  return
 else:
  time=o.realtime
  porosity=voxelPorosity(resolution=600,start=aabbExtrema()[0],end=aabbExtrema()[1])
  Maxoverlap=utils.maxOverlapRatio()
  stress=getStress().trace()/3.
  row_contents = [time,time,porosity,Maxoverlap,stress]
  append_list_as_row('PH102_Haustein_data.csv', row_contents)

Hi,

Thank you for your reply.
Sorry, I have attached the whole code. I want to have to contact laws, one for the contact between particles and the walls and the other on for the contact between the particles.
The question is how can add to contact law in one engine and specify which law is used for particle-particle interaction and which for particle-wall interaction?

> I want to have to contact laws

to = two?

As I mentioned before, in principle it is possible to have multiple Law2 functors. Just specify them the same way as you have multiple Ig2 functors.

The other point is how it will work.
Law2 to be used for the interaction is determined by IGeom and IPhys.

IGeom is created / updated by Ig2 functor.
Ig2 to be used is determined by shapes of interacting particles.
E.g., the collider detects interaction between sphere-sphere, so Ig2_Sphere_Sphere_ScGeom is used, resulting in ScGeom
Or e.g. the interaction is between sphere-facet, so Ig2_Facet_Sphere_ScGeom is used, resulting again in ScGeom (but can be other IGeoms in principle).

Similarly it works for IPhys
IGeom is created / updated by Ip2 functor.
Ip2 to be used is determined by materials of interacting particles.
E.g. the materials of both interaction bodies is LudingMat, then you can use Ip2_LudingMat_LudingMat_LudingPhys to create LudingPhys and together with ScGeom you can use Law2_ScGeom_LudingPhys_Basic.
Or e.g. the materials of both interaction bodies is FrictMat, then you can use Ip2_FrictMat_FrictMat_MindlinPhys to create MindlinPhys and together with ScGeom you can use Law2_ScGeom_MindlinPhys_Mindlin.

Now, you could use LudingMat for spheres and FrictMat for walls.
Then you could use Ip2_LudingMat_LudingMat_LudingPhys and Ip2_FrictMat_LudingMat_MidlinPhys.
Sphere-sphere (or more precisely Luding-Luding) contanct would result in LudingPhys and Law2_ScGeom_LudingPhys_Basic would be used.
Sphere-wall (or more precisely FrictMat-LudingMat) contacts would result in MindlinPhys and Law2_ScGeom_MindlinPhys_Mindlin would be used.
However, Ip2_FrictMat_LudingMat_MidlinPhys is not implemented in Yade :-(
You would have to implement it in the C++ source code.

How it works for the Cpm case I described before:
The point here is that CpmMat is derived from FrictMat, so CpmMat may be treated as FrictMat.
Spheres have CpmMat, walls have FrictMat.
Ip2 are Ip2_CpmMat_CpmMat_CpmPhys and Ip2_FrictMat_FrictMat_FrictPhys
Sphere-sphere (or more precisely cpm-cpm) contanct would result in CpmPhys and Law2_ScGeom_CpmPhys_Cpm would be used.
Sphere-wall (or more precisely cpm-frict wich are essentially frict-frict) contacts would result in FrictPhys and Law2_ScGeom_FrictPhys_CundallStrack would be used.

You cannot do the same trick with LudingMat, as LudingMat is derived directly from Material class, not FrictMat... [3]

A general note aboutsing multiple shapes / materials.
You have to be careful:
1) not to miss a desired combination.
E.g. without Ig2_Facet_Sphere_ScGeom, spheres would pass through facets without interaction, but without any error at the same time
2) not to define ambiguous functors - the same input could result in more than one result. Yade woudl give you run time error in that case

cheers
Jan

[3] https://yade-dem.org/doc/yade.wrapper.html#material