yade-users team mailing list archive

[Jan Stránský <question252656@xxxxxxxxxxxxxxxxxxxxx>]

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

Jan Stránský proposed the following answer:
Hi Colin,
sorry for late reply. I tried your script, but to run it, I also
need ccore_yade module. Could you please send it in email or place it
somewhere to the internet and send a link?
thanks
Jan



2014-08-05 16:51 GMT+02:00 Colin Power <question252656@xxxxxxxxxxxxxxxxxxxxx
>:

> Question #252656 on Yade changed:
> https://answers.launchpad.net/yade/+question/252656
>
> Colin Power posted a new comment:
> Yes, there are multiple facets in the simulation and there are also are
> multiple material types. I am also using a buoyancy model that is
> outlined in this example:
>
> https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps
> /apply-buoyancy-clumps.py
>
> I have pasted my script below. The ccore_yade module just contains
> helper functions for file management.
>
> import sys, time, random, os, gts, math
> from yade import ymport
> from yade import pack
> from yade import qt
> import ccore_yade
>
> # Simulation Parameters
> nbIter=6000000
> filename = 'seabed'
> volumeName = 'keel'
> iterPeriod = 100
>
> folder = ccore_yade.newSimulationFolder(filename)
> print "Simulation is being stored in: " + folder
> filepath = folder+'/'+filename
>
> #used for randomDensePack, experimenting with this.
> memoizeDb='keel_packings.sqlite'
>
> #define material properties:
> shearModulus                    = 3.6e9
> poissonRatio                    = 0.28
> youngModulus                    = 2*shearModulus*(1+poissonRatio)
> angle                                   = radians(36) #friction angle in
> radians
> angle_r                                 = radians(36) #friction angle in
> radians
> rho_p                                   = 917         #density of particles
> rho_f                                   = 1029        #density of the
> fluid, rho_f > rho_p = floating particles
> rho_r                   = 1602        #density of seafloor rock
> rho_s                   = 8050        #density of steel
> v_iw                                    = 0           #velocity of
> increasing water-level
> gravity                 = -9.80665    #m/s
> damping                 = 0.3         #ratio
> radius                  = 0.05        #m
> seabed_velocity         = 0.05        #m/s
>
> #define water boundaries:
> boundaryMin = Vector3(-6, -6, -0.14)
> boundaryMax = Vector3(10, 10, 10)
> waterLevel = boundaryMin[2]
> t0 = O.time
>
> #define colors:
> sphereColor = (0.0,0.0,0.05)  #dirty yellow
> waterColor      = (.2,.2,.7)  #blue
>
> #define materials:
> id_Mat=O.materials.append(CpmMat(young=youngModulus, poisson=poissonRatio,
> density=rho_p, frictionAngle=angle, sigmaT=5.5e3, relDuctility=30,
> epsCrackOnset=1e10, isoPrestress=0, damLaw=0) )
> Mat=O.materials[id_Mat]
>
> id_Mat=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_r,
> frictionAngle=angle_r))
> rock=O.materials[id_Mat]
>
> #steel materials have been seperated for filtering in Paraview
> steel_Mat1=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel1=O.materials[steel_Mat1]
>
> steel_Mat2=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel2=O.materials[steel_Mat2]
>
> steel_Mat3=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel3=O.materials[steel_Mat3]
>
> # Import mesh files to create facets
> seabed = O.bodies.append(ymport.gmsh(meshfile=filename+'.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.75,0.35,0), material=rock, wire=False))
> leftWall = O.bodies.append(ymport.gmsh(meshfile='leftWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel1, wire=False))
> rightWall = O.bodies.append(ymport.gmsh(meshfile='rightWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel2, wire=False))
> topWall = O.bodies.append(ymport.gmsh(meshfile='topWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel3, wire=False))
>
> # Create water level indicator
> idsWaterFacets =  []
> idsWaterFacets.append(O.bodies.append(facet( \
>                         [ boundaryMin, [boundaryMax[0], boundaryMin[1],
> boundaryMin[2]], [boundaryMax[0], boundaryMax[1], boundaryMin[2]] ], \
>                         fixed=True, material=FrictMat(young=0),
> color=waterColor, wire=False))) #no interactions will appear
> idsWaterFacets.append(O.bodies.append(facet( \
>                         [ [boundaryMax[0], boundaryMax[1],
> boundaryMin[2]], [boundaryMin[0], boundaryMax[1], boundaryMin[2]],
> boundaryMin ], \
>                         fixed=True, material=FrictMat(young=0),
> color=waterColor, wire=False))) #no interactions will appear
>
> surface=gts.read(open(volumeName+'.gts'))
> volume=pack.inGtsSurface(surface)
>
> O.bodies.append(pack.randomDensePack(predicate=volume, radius=radius,
> material=Mat, cropLayers=0, rRelFuzz=0.6666, spheresInCell=50000,
> memoizeDb=None, useOBB=False, memoDbg=False, color=None,
> returnSpherePack=False))
>
> #Define timestep
> O.dt=0.5*PWaveTimeStep()
>
> ccore_yade.recordTimeStep(filepath, iterPeriod)
>
> O.engines=[
>         ForceResetter(),
>
>         InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
>         InteractionLoop(
>                 [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
>                 [Ip2_Mat_CpmMat_CpmPhys(cohesiveThresholdIter=10000),
> Ip2_FrictMat_CpmMat_FrictPhys()],
>                 [Law2_ScGeom_FrictPhys_CundallStrack(),
> Law2_ScGeom_CpmPhys_Cpm(omegaThreshold=1, yieldEllipseShift=0,
> yieldSurfType=0)]
>         ),
>
>         NewtonIntegrator(damping=damping, gravity=[0, 0, gravity],
> label='integrator'),
>     TranslationEngine(dead=False, ids=seabed, translationAxis=[-1,0,0],
> velocity=seabed_velocity, label='seabed_movement'),
>
>     #Moveable walls for future simulations
>     #TranslationEngine(dead=False, ids=leftWall, translationAxis=[0,1,0],
> velocity=0.1, label='leftWall_movement'),
>     #TranslationEngine(dead=False, ids=rightWall,
> translationAxis=[0,-1,0], velocity=0.1, label='rightWall_movement'),
>     #TranslationEngine(dead=False, ids=topWall, translationAxis=[0,0,-1],
> velocity=0.1, label='topWall_movement'),
>
>     # Save data for Paraview
>     VTKRecorder(fileName=filepath+'_vtk', recorders=['all', 'cpm'],
> iterPeriod=iterPeriod),
>
>     PyRunner(command='O.pause()', iterPeriod=nbIter),
>     PyRunner(iterPeriod=25, command='applyBuoyancy()', label='buoLabel'),
>     PyRunner(iterPeriod=10000, command='print O.iter')
> ]
>
> print "Start simulation: " + filename
> qt.View() # launch Yade with the 3D view
> O.pause() # launch Yade with the simulation paused
>
> # Based on
> https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps/apply-buoyancy-clumps.py
> def applyBuoyancy():
>         global waterLevel
>         waterLevel = (O.time-t0)*v_iw + boundaryMin[2]
>         for b in O.bodies:
>                 zMin = 1e9
>                 zMax = -1e9
>                 if b.isStandalone and isinstance(b.shape,Sphere):
>                         rad = b.shape.radius
>                         zMin = b.state.pos[2] - rad
>                         dh = min((waterLevel - zMin), 2*rad)    #to get
> sure, that dh is not bigger than 2*radius
>                 elif b.isClump:                         #determine rad,
> zMin and zMax for clumps:
>                         for ii in b.shape.members.keys():
>                                 pos = O.bodies[ii].state.pos
>                                 zMin =
> min(zMin,pos[2]-O.bodies[ii].shape.radius)
>                                 zMax =
> max(zMax,pos[2]+O.bodies[ii].shape.radius)
>                         #get equivalent radius from clump mass:
>                         rad = (
> 3*b.state.mass/(4*pi*O.bodies[b.shape.members.keys()[0]].mat.density)
> )**(1./3.)
>                         #get dh relative to equivalent sphere, but acting
> when waterLevel is between clumps z-dimensions zMin and zMax:
>                         dh = min((waterLevel - zMin)*2*rad/(zMax - zMin),
> 2*rad)
>                 else:
>                         continue
>                 if dh > 0:
>                         F_buo = -1*(pi/3)*dh*dh*(3*rad -
> dh)*rho_f*integrator.gravity
>                            # = -V*rho*g
>                         O.forces.addF(b.id, F_buo, permanent=True)
>
> --
> 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.