yade-users team mailing list archive

[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.