yade-users team mailing list archive

[Christian Jakob <question230621@xxxxxxxxxxxxxxxxxxxxx>]

New question #230621 on Yade:
https://answers.launchpad.net/yade/+question/230621

Hi,

I want to model increasing water level in a sandy soil. Therefore I want to include buoyancy forces, that should act on particles (not at contacts).

I tried O.forces.addF() method, but it just adds forces for one single step. So I tried ForceEngine and for one single particle, the script is working fine:

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

#define model properties:
shear_modulus			= 3.2e10
poisson_ratio			= 0.15
young_modulus			= 2*shear_modulus*(1+poisson_ratio)
friction_coeff			= 10
angle					= atan(friction_coeff)
rho_p					= 2650	#density of particles

v_iw					= 1  #velocity of increasing water-level

#material:
id_Mat=O.materials.append(FrictMat(young=young_modulus,poisson=poisson_ratio,density=rho_p,frictionAngle=angle))
Mat=O.materials[id_Mat]

#create sphere and box:
O.bodies.append(sphere(center=(0,0,1), radius=1, fixed=True, material=Mat))
O.bodies.append(box(center=(0,0,-.1), extents=(1,1,.1), fixed=True, material=Mat))

#define engines:
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.0,betas=0.0,label='ContactModel')],
		[Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False,label='Mindlin')]
		),
	NewtonIntegrator(gravity=(0,0,0),damping=0.0,label='integrator')
]
#initiate force engines for all particles:
for b in O.bodies:
	if (b.isStandalone and isinstance(b.shape,Sphere)) or b.isClump:
		O.engines=O.engines+[ForceEngine(ids=[b.id],label='fe%i'%b.id)]

def apply_buo(water_height,saturatedList):
	for b in O.bodies:
		if b not in saturatedList:
			if b.isStandalone and isinstance(b.shape,Sphere):
				pos = b.state.pos
				rad = b.shape.radius
				Id = b.id
			else:
				continue
			h = pos[2] - rad
			if h < water_height:
				dh = water_height - h
				dh = min(dh,2*rad)	#to get sure, that dh is not bigger than 2*radius
				F_buo = Vector3(0,0,(pi/3)*dh*dh*(3*rad - dh)*9810)	#=V*rho*g
				
				################################# TRICKY PART HERE 
				#O.forces.addF(Id,F_buo)	#not working, it remembers the force just for 1 step
				fe0.force = F_buo			#can not be used for many particles ... hmpf

				if (h+2*rad) < water_height:
					saturatedList.append(O.bodies[Id])
					print 'particle is fully surrounded by water'
	return saturatedList

saturatedList = []
t0 = O.time
O.dt=1e-4

O.engines=O.engines+[PyRunner(iterPeriod=10,command='saturatedList=apply_buo(((O.time-t0) * v_iw),saturatedList)',label='buoLabel')]
O.engines=O.engines+[PyRunner(iterPeriod=999,command='print O.forces.f(0)',label='fLabel')]

O.run(25000,True)

#theoretical buoyancy force:
print 'theoretical F_buo = ',4*pi*9810/3

###endofscript

The problem occurs, when I have thousands of particles in the model. How can I find out which particle has which ForceEngine? Do you have other ideas how to solve this problem?

Regards,

Christian

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