yade-users team mailing list archive

[Anton Gladky <question230621@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Anton Gladky proposed the following answer:
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 <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.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp

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