yade-users team mailing list archive

[Christian Jakob <question232941@xxxxxxxxxxxxxxxxxxxxx>]

Question #232941 on Yade changed:
https://answers.launchpad.net/yade/+question/232941

Christian Jakob proposed the following answer:
You can also investigate the script I wrote for getting a force-distance-plot.
There data is stored in python lists and plotted with pylab/matplotlib.
If you want, you can modify this script and add the displacement O.bodies[X].state.displ(), like Jan suggested.

hope it helps,

christian

#!/usr/bin/python
# -*- coding: utf-8 -*-

### properties:
density			= 2650
shear_modulus 	= 2e5
poisson_ratio 	= 0.15
young_modulus	= 2*shear_modulus*(1+poisson_ratio)
suction 		= 3000
velocity		= 0.1	#in x-direction
a=True

### define a material:
mat = FrictMat(young=young_modulus,poisson=poisson_ratio,density=density,frictionAngle=1)
O.materials.append(mat)

### append 2 spheres:
R=.0003
d=.000015	#positive for overlaps, negative for distant
O.bodies.append(sphere([0,0,0],radius=R, material = mat, fixed=True))#id 0
O.bodies.append(sphere([2*R-d,0,0],radius=R, material = mat, fixed=True))#id 1

### define engines:
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_MindlinCapillaryPhys(betan=0,betas=0,label='ContactModel')],
		[Law2_ScGeom_MindlinPhys_Mindlin(neverErase=True,label='Mindlin')]
		),
	Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=suction,label='CapPhys'),
	NewtonIntegrator(damping=0,label='integrator'),
	PyRunner(iterPeriod=1,command='a=getData(a)')
]
#CapPhys.dead=True
O.bodies[1].state.vel=Vector3(velocity,0,0)
O.dt=1e-6

capForce=[]
normForce=[]
normCapForce=[]
distance=[]

def getData(a):
	distance.append(O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0]-2*O.bodies[0].shape.radius)
	try:
		i = O.interactions[0,1]
		if i.isReal:
			capForce.append(-1*i.phys.fCap[0])
			normForce.append(i.phys.normalForce[0])
			normCapForce.append(i.phys.normalForce[0]-i.phys.fCap[0])
		else:
			if a:
				print 'interaction lost at step',O.iter
				a=False
			capForce.append(0)
			normForce.append(0)
			normCapForce.append(0)
	except IndexError:
		capForce.append(0)
		normForce.append(0)
		normCapForce.append(0)
	return a

if d <= 0.0:
	CapPhys.createDistantMeniscii=True
	O.run(1,True)
	CapPhys.createDistantMeniscii=False

O.run(600,True)


### make nice figure:

from pylab import *		#overwrites angle, box, ... see 0-generate.py
rc('text',usetex=True)
rc('font',**{'family':'serif','serif':['Computer Modern']})

from matplotlib.font_manager import FontProperties
font = FontProperties()
font.set_size('x-large')

figure()
axhline(color='gray')
axvline(color='gray')
plot(distance, normCapForce, label = 'normal + capillary force')
plot(distance, normForce, 'g--', label = 'normal force')
plot(distance, capForce, 'r--', label = 'capillary force')
xlabel('distance in [$m$]',fontproperties=font)
ylabel('force in [$N$]',fontproperties=font)
ylim([-0.00015,0.00025])
xlim([-0.00002,0.00005])
#grid()
legend()
#savefig('force-dist-plot-cap+norm.png')
show()

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.