# yade-users team mailing list archive

## Re: [Question #696052]: Two contact law

 Thread Previous • Date Previous • Date Next • Thread Next

```Question #696052 on Yade changed:

Mithushan Soundaranathan posted a new comment:
#!/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)

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

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

######################################################################

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

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)

--