apply buoyancy forces on particles

Hi,

from my point of view it is a bad idea to use just 1 engine for 1
particle. You can create a forceBuoyancyEngine, which will
accept the list of affected particles and calculate forces. Please,
include also in the code references on equation in the Paper.

Cheers,

Anton

2013/6/12 Christian Jakob <email address hidden>:
> 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.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Thanks anton, i will think about that.

Why not this:

if h < water_height:
   b.state.mass = 4*pi/3*rad**3*rho_s - pi/3*dh*dh*(3*rad - dh)*rho #so that apparent weight will be solid weight - buoyancy

Wow, a very nice idea. Easy to implement and fast in calculation ... *thumbs up* ^^

Thank you very much!

Hi, Bruno
Is there sometime wrong with your sript? When we calculate inertia force, the rho_s is used. but you change the b.state.mass , is the (rho_s-rho_fluid) is used?
Maybe I am wrong because I donnot know what is the dh means?

Hi,Anton
Using just 1 engine for 1particle may cause the script too slow, but it can mantain the force for some iterPeriod for different particle without addForce every iter. What is more effective, I donnot know.

        ForceEngine(ids=ids_spSlurry,force=(0,0,mass_p*acc_grav*density_fluid/density_p),label="buoyForce"),

I must admit my trick has the drawback of changing the dynamic of the system:
in acceleration=force/mass, mass will be changed while it should not (only force should be different)

If necessary, it can be fixed using b.state.densityScaling (see Newton.cpp:179):
b.state.mass = apparentMass
b.state.densityScaling = realMass/apparentMass.

And: newton.densityScaling=True

In that case acceleration is correct again.

Liucheng, I think your question is for Christian, not for me. I simply duplicated his formula in a different form.
I guess you assume each particle is surrounded by fluid, Christian does not.

>If necessary, it can be fixed using b.state.densityScaling (see Newton.cpp:179):
>b.state.mass = apparentMass
>b.state.densityScaling = realMass/apparentMass.

ok, thanks again, i will fix that.

btw, is it worth to create a small script about buoyancy in examples section?

2013/6/14 Christian Jakob <email address hidden>:
>
> btw, is it worth to create a small script about buoyancy in examples
> section?

Sure.

Anton

> b.state.densityScaling = realMass/apparentMass

if i decrease mass m, then accel. a = F/m should increase...
but densityScaling increases accel. again:

if (densityScaling) linAccel*=state->densityScaling;

so it should be

b.state.densityScaling = apparentMass/realMass

right?

You are right, sorry :)