yade-users team mailing list archive

[Daria <question692907@xxxxxxxxxxxxxxxxxxxxx>]

New question #692907 on Yade:
https://answers.launchpad.net/yade/+question/692907

I have got case in which flow around cylinder is simulated. This is YADE-OpenFOAM coupling case. I have got closed cylinder, but at some time particles enter it and simulation fails. How can I fix it?

scriptYade.py
#####
from __future__ import print_function
import sys
from yadeimport import *
from yade.wrapper import *
from yade.utils import *

initMPI()                           #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling();     #Initialize the engine
fluidCoupling.getRank();            #part of Initialization.


#example of spheres in shear flow : two-way point force coupling
class simulation():

	def __init__(self):
		O.periodic = True
		O.cell.setBox(1,0.4,0.2)

		numspheres=1000
		young = 1e10#5e6
		young_steel = 2e11
		density = 2#2000
		density_steel = 7700
		poisson = 0.3#0.5
		normalCohesion = 1e7
		shearCohesion = 1e7
		frictionAngle = radians(15)

		mat1 = CohFrictMat(normalCohesion=normalCohesion, shearCohesion=shearCohesion, young=young, poisson=poisson, frictionAngle=frictionAngle, density=density, momentRotationLaw=True, label='spheremat')
		O.materials.append(mat1)
		mat2 = CohFrictMat(normalCohesion=normalCohesion, shearCohesion=shearCohesion, young=young_steel, poisson=poisson, frictionAngle=frictionAngle, density=density_steel, momentRotationLaw=True, label='wallmat')
		O.materials.append(mat2)

		epsilon = 1e-08
		minval = 0 + epsilon
		maxx = 1.0 - epsilon
		maxy = 0.4 - epsilon
		maxz = 0.2 - epsilon
		#wall coords, use facets for wall BC:
		v0 = Vector3(minval, minval, minval)
		v1 = Vector3(minval,minval,maxz)
		v2 = Vector3(maxx,minval,minval)
		v3 = Vector3(maxx,minval,maxz)

		v4 = Vector3(minval,maxy,minval)
		v5 = Vector3(minval,maxy,maxz)
		v6 = Vector3(maxx,maxy,minval)
		v7 = Vector3(maxx, maxy, maxz)

		lf0 = facet(vertices=[v0,v1,v2], material='wallmat')
		O.bodies.append(lf0)
		lf1 = facet(vertices=[v1,v2,v3], material='wallmat')
		O.bodies.append(lf1)

		uf0 = facet(vertices=[v4,v5,v6], material='wallmat')
		O.bodies.append(uf0)
		uf1 = facet(vertices=[v5,v6,v7], material='wallmat')
		O.bodies.append(uf1)

		ff0 = facet(vertices=[v0,v2,v6], material='wallmat')
		O.bodies.append(ff0)
		ff1 = facet(vertices=[v0,v6,v4], material='wallmat')
		O.bodies.append(ff1)

		bf0 = facet(vertices=[v1,v3,v7], material='wallmat')
		O.bodies.append(bf0)
		bf1 = facet(vertices=[v1,v7,v5], material='wallmat')
		O.bodies.append(bf1)

		oriBody = Quaternion(Vector3(1,0,0),(0))
		radius = 0.05
		O.bodies.append(geom.facetCylinder((0.5,0.2,0.1),radius=radius,height=0.2,orientation=oriBody,segmentsNumber=100,wallMask=4,material='wallmat'))

		#spheres
		mn, mx= Vector3(minval + epsilon, minval + epsilon, minval + epsilon), Vector3(maxx / 2 - radius - epsilon, maxy - epsilon, maxz - epsilon)
		sp = pack.SpherePack();
		sp.makeCloud(mn,mx,rMean=0.003,rRelFuzz=0.5, num=numspheres)
		O.bodies.append([sphere(center,rad,material='spheremat') for center,rad in sp])
		mn, mx= Vector3(maxx / 2 + radius + epsilon, minval + epsilon, minval + epsilon), Vector3(maxx - epsilon, maxy - epsilon, maxz - epsilon)
		sp = pack.SpherePack();
		sp.makeCloud(mn,mx,rMean=0.003,rRelFuzz=0.5, num=numspheres)
		O.bodies.append([sphere(center,rad,material='spheremat') for center,rad in sp])
		sphereIDs = [b.id for b in O.bodies if type(b.shape)==Sphere]

		#coupling engine settings
		fluidCoupling.setNumParticles(len(sphereIDs))
		fluidCoupling.setIdList(sphereIDs)
		fluidCoupling.isGaussianInterp=False;  #use pimpleFoamYade for gaussianInterp

		# Integrator
		newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
		# add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.

		O.timingEnabled==True

		O.engines=[
			ForceResetter(),
			InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
			InteractionLoop(
				[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
				[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True)],
				[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()],
			GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.3, label = "ts"),
			fluidCoupling, #to be called after timestepper
			PyRunner(command='sim.printMessage()', iterPeriod= 1, label='outputMessage'),
			VTKRecorder(fileName='yadep/3d-vtk-',recorders=['spheres', 'facets', 'intr', 'force'],virtPeriod=0.001)
		]

	def printMessage(self):
		print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
		print("********************************YADE-TIME = " + str(O.time) +" **********************************")

	def printState(self, n):
		print('blaaaaa ' + str(n))

	def printDt(self):
		print(O.dt)

	def setDt(self, dt):
		O.dt = dt

	def irun(self,num):
		O.run(num,1)

if __name__=="__main__":
	sim = simulation()
	sim.irun(10000000)
	fluidCoupling.killMPI()

import builtins
builtins.sim=sim

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.