cpm_crash_test.py

import sys, time, random, os, gts, math
from yade import ymport
from yade import pack
from yade import qt

#Simulation Parameters
nbIter=100000000000
iterPeriod = 100
filename = "crash_test"

#define material properties:
shearModulus            = 3.6e9
poissonRatio            = 0.28
youngModulus            = 2*shearModulus*(1+poissonRatio)
angle                   = 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
v_iw                    = 0           #velocity of increasing water-level
gravity                 = -9.80665    #m/s
damping                 = 0.3         #ratio

#define water boundaries:
boundaryMin = Vector3(-10, -10, -0.14)
boundaryMax = Vector3(10, 10, 10)
waterLevel = boundaryMin[2]
t0 = O.time

#define colors:
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]

# 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('keel.gts'))
volume=pack.inGtsSurface(surface)

O.bodies.append(pack.randomDensePack(predicate=volume, radius=0.05, material=Mat, cropLayers=0, rRelFuzz=0.6666, spheresInCell=5000, memoizeDb=None, useOBB=False, memoDbg=False, color=None, returnSpherePack=False))


#Define timestep
O.dt=0.5*PWaveTimeStep()

O.engines=[
    ForceResetter(),

    InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_CpmMat_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'),

    # 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)