Question #680892 on Yade changed:
https://answers.launchpad.net/yade/+question/680892
Status: Answered => Open
Francesco is still having a problem:
Hi,
I attached the code. I commented some functions related to data storing
that could create problems.
Thank you so much.
Francesco
Code:
from numpy import linspace
from yade import pack, qt
import numpy as np
import os
import glob
import re
import smtplib
from email.MIMEMultipart import MIMEMultipart
from email.MIMEText import MIMEText
from email.MIMEBase import MIMEBase
from email import encoders
######################################################
# Parameters
######################################################
mi=10**-3 # Millimetres transformation
mu=10**-6 # Micron transformation
x1=10*mi # Box half side
z1=10*mi # Box heigth
dist=100*mi # Nozzle target distance
h=z1+(dist) # Nozzle heigth from z=0
rIn=6*mi # Inner radius
rOut=12*mi # Cylinder radius
eps=0.1*mi # Tollerance
itereps=5
rPack=rOut-eps # Pack radius
rd=rOut-rIn # Difference between radius
l2=80*mi # Cylinder length
l_cloud=50*mi # Particles cloud length
c1=(4*l_cloud)/9 # External particles pack range
c2=(1*l_cloud)/9 # Middle particles pack range
l3=1*mi # Distance from the nozle indentation
v=57 # Particles initial velocity
vi_y=-v*cos(pi*angle/180)
vi_z=-v*sin(pi*angle/180)
vf1=57 # Fluid velocity 1
vf2=220 # Fluid velocity 2
rMean=40*mu # Mean particles radius
Rho=3.61 # Air density 1
Cd=0.47 # Sphere drag coefficient
wt=100000 # Box rotation
angle=22.5 # Nozzle angle
box_angle=45
vx=10000
vz=10000
######################################################
# Nozzle creation
######################################################
delta=((pi*0.5)/16)
poly=((rIn,h),
(rOut-(rd*cos(delta)),h+(rIn*sin(delta))),
(rOut-(rd*cos(2*delta)),h+(rIn*sin(2*delta))),
(rOut-(rd*cos(3*delta)),h+(rIn*sin(3*delta))),
(rOut-(rd*cos(4*delta)),h+(rIn*sin(4*delta))),
(rOut-(rd*cos(5*delta)),h+(rIn*sin(5*delta))),
(rOut-(rd*cos(6*delta)),h+(rIn*sin(6*delta))),
(rOut-(rd*cos(7*delta)),h+(rIn*sin(7*delta))),
(rOut-(rd*cos(8*delta)),h+(rIn*sin(8*delta))),
(rOut-(rd*cos(9*delta)),h+(rIn*sin(9*delta))),
(rOut-(rd*cos(10*delta)),h+(rIn*sin(10*delta))),
(rOut-(rd*cos(11*delta)),h+(rIn*sin(11*delta))),
(rOut-(rd*cos(12*delta)),h+(rIn*sin(12*delta))),
(rOut-(rd*cos(13*delta)),h+(rIn*sin(13*delta))),
(rOut-(rd*cos(14*delta)),h+(rIn*sin(14*delta))),
(rOut-(rd*cos(15*delta)),h+(rIn*sin(15*delta))),
(rOut,h+rIn),(rOut,h+rIn+(l2/8)),(rOut,h+rIn+(l2/4)),
(rOut,h+rIn+((3*l2)/8)),(rOut,h+rIn+((4*l2)/8)),(rOut,h+rIn+((5*l2)/8)),
(rOut,h+rIn+((6*l2)/8)),(rOut,h+rIn+((7*l2)/8)),(rOut,h+rIn+(l2))
)
thetas=linspace(0,2*pi,num=32,endpoint=True)
Nozzle=O.materials.append(FrictMat(density=7.8*10**3,young=210e9,poisson=0.3,label="Nozzle"))
pts=pack.revolutionSurfaceMeridians([[(pt[0],pt[1]) for pt in poly] for theta
in thetas],thetas)
surf=pack.sweptPolylines2gtsSurface(pts)
nozzle=O.bodies.append(pack.gtsSurface2Facets(surf))
######################################################
# Creating Target
######################################################
Aluminum=O.materials.append(FrictMat(density=2.81*10**3,young=71e9,poisson=0.33,label="Aluminum"))
p=utils.box((x1,0,0),(x1,x1,z1),fixed=True,material="Aluminum")
box=O.bodies.append(p)
l=len(O.bodies)
t_id=l-1
p1=utils.box((-x1,0,0),(x1,x1,z1),fixed=True,material="Aluminum")
box1=O.bodies.append(p1)
l_1=len(O.bodies)
t_id1=l_1-1
p2=utils.box((2*x1,0,0),(2*x1,x1,z1),fixed=True,material="Aluminum")
box2=O.bodies.append(p2)
l_2=len(O.bodies)
t_id2=l_2-1
p3=utils.box((-2*x1,0,0),(2*x1,x1,z1),fixed=True,material="Aluminum")
box3=O.bodies.append(p3)
l_3=len(O.bodies)
t_id3=l_3-1
p4=utils.box((x1,-2*x1,0),(x1,x1,z1),fixed=True,material="Aluminum")
box4=O.bodies.append(p4)
l_4=len(O.bodies)
t_id4=l_4-1
######################################################
# Creating shots
######################################################
V_cloud=((2*rOut)*(2*rOut))*l_cloud*1.1
V_cyl=(pi*rOut*rOut*(l_cloud))
n_shot=1287
dens=(n_shot/V_cyl)
n_shot_cloud=int(dens*V_cloud)
corner=(rOut)
Ceramic=O.materials.append(FrictMat(density=2.3*10**3,young=300e9,poisson=0.30,label="Ceramic"))
sp=pack.SpherePack()
sp.makeCloud((-corner,-corner,h+rIn+l3),(corner,corner,h+rIn+l3+(1.1*l_cloud)),num=n_shot_cloud,psdSizes=[53*10**-6,63*10**-6,90*10**-6,106*10**-6,125*10**-6],psdCumm=[0,0.42,0.76,0.92,1])
pred_1=pack.inCylinder((0,0,h+rIn+l3),(0,0,h+rIn+l3+(c1)),radius=rPack)
pred_2=pack.inCylinder((0,0,h+rIn+l3+c1),(0,0,h+rIn+l3+(c1+c2)),radius=rPack)
pred_3=pack.inCylinder((0,0,h+rIn+l3+c1+c2),(0,0,h+rIn+l3+(c1+c2+c1)),radius=rPack)
sp2=pack.filterSpherePack(pred_1,sp,returnSpherePack=True)
sp3= pack.filterSpherePack(pred_2,sp,returnSpherePack=True)
sp4=pack.filterSpherePack(pred_3,sp,returnSpherePack=True)
sp2.toSimulation()
j1=len(O.bodies)
sp3.toSimulation()
j2=len(O.bodies)
sp4.toSimulation()
j=len(O.bodies)
pack_rotation=range(l_4+1,j)
O.dt=utils.PWaveTimeStep()
([utils.setBodyVelocity(i,(0,vi_y,vi_z),axis="xyz") for i in
range(l_4,j)])
#periodrecord=20
#dzrecord=1*mi
#nrecord=int((l_cloud+rIn+l3)/dzrecord)
#startrecordtime=0
#endrecordtime=int((1*mi)/(O.dt*v))
#totalrecordtime=endrecordtime*nrecord
#deltarecordtime=(endrecordtime-startrecordtime)
#ndorecord=(int(deltarecordtime/periodrecord))*nrecord*2
#deltarecord=(endrecordtime-startrecordtime)*O.dt
#firstrecord=40
#lastrecord=45
#hin=h+rIn+l3
#shot=[]
#mass=[]
#Vnorm=[]
#def txt():
# for ii in range(l_4,j):
# if O.bodies[ii].state.pos[2]>h:
# if O.bodies[ii].state.pos[2]<hin:
# vb=O.bodies[ii].state.vel[2]
# radius=O.bodies[ii].shape.radius
# Fd=0.5*Rho*Cd*(pi*(radius)**2)*(vb-(-vf2))**2
# O.forces.setPermF(ii,(0,0,-Fd))
# else:
# if O.bodies[ii].state.vel[2]>0:
# vb=O.bodies[ii].state.vel[2]
# radius=O.bodies[ii].shape.radius
#
Fd=0.5*Rho*Cd*(pi*(radius)**2)*(vb-(-vf2))**2
# O.forces.setPermF(ii,(0,0,-Fd))
# else:
# O.forces.setPermF(ii,(0,0,0))
# if O.bodies[ii].state.pos[2] < h:
# if O.bodies[ii].state.pos[2] > h-eps:
#
name="/home/michelangelo/Scrivania/YADE_PIONA/NOZZLE_FLOW_EDGE_37/RESULTS_2/velocity"+str(ii)+".txt"
# with open(name,"a") as f:
# if os.stat(name).st_size == 0 :
# vout=O.bodies[ii].state.vel
#
V=(((vout[0])**2)+((vout[1])**2)+((vout[2])**2))**(0.5)
# z=O.bodies[ii].state.pos[2]
# m=O.bodies[ii].state.mass
#
f.write("{:.2f}\t{:.2f}\t{:.2f}\t{:.3f}\t{}\n".format(vout[0],vout[1],vout[2],z,ii))
# mass.append(m)
# Vnorm.append(V)
# if
O.time>(firstrecord*deltarecord):
# if
O.time<(lastrecord*deltarecord):
# shot.append(ii)
# else:
# pass
# else:
# pass
# else:
# pass
#start=int((dist)/(O.dt*v))
#
#rmax=((125*10**-6)/2)
#
#def txt_3():
# global yref0
# global vref0
# for ii in shot:
#
name="/home/michelangelo/Scrivania/YADE_PIONA/NOZZLE_FLOW_EDGE_37/RESULTS_3/contact"+str(ii)+".txt"
# with open(name,"a+") as f:
# if O.bodies[ii].state.vel[2]<-1:
# if os.stat(name).st_size == 0 :
# if O.bodies[ii].state.pos[0]<0:
# if
O.bodies[ii].state.pos[2]<z1+(O.bodies[ii].state.pos[0]+rmax)*tan((pi*(angle))/180):
#
cp=O.bodies[ii].state.pos
#
vc=O.bodies[ii].state.vel
# t=O.time
#
f.write("{:.6f}\t{:.6f}\t{:.6f}\t{:.2f}\t{:.2f}\t{:.2f}\t{}\t{}\n".format(cp[0],cp[1],cp[2],vc[0],vc[1],vc[2],ii,t))
# else:
# pass
# else:
# if
O.bodies[ii].state.pos[2]<(z1+rmax)-(O.bodies[ii].state.pos[0])*tan((pi*(90-angle))/180):
#
cp=O.bodies[ii].state.pos
#
vc=O.bodies[ii].state.vel
# t=O.time
#
f.write("{:.6f}\t{:.6f}\t{:.6f}\t{:.2f}\t{:.2f}\t{:.2f}\t{}\t{}\n".format(cp[0],cp[1],cp[2],vc[0],vc[1],vc[2],ii,t))
# else:
# pass
# else:
# if O.bodies[ii].state.pos[0]<0:
# if
O.bodies[ii].state.pos[2]<z1+(O.bodies[ii].state.pos[0]+rmax)*tan((pi*(angle))/180):
# lines=f.readlines()
# for x in lines:
#
old=(re.split('\t|\n',x))
#
yref0=float(old[1])
#
vref0=float(old[4])
# if
yref0-mi<=round(O.bodies[ii].state.pos[1],6)<=yref0+mi:
# pass
# else:
# if
round(O.bodies[ii].state.vel[1],2)==vref0:
# pass
# else:
#
cp=O.bodies[ii].state.pos
#
vc=O.bodies[ii].state.vel
# t=O.time
#
f.write("{:.6f}\t{:.6f}\t{:.6f}\t{:.2f}\t{:.2f}\t{:.2f}\t{}\t{}\n".format(cp[0],cp[1],cp[2],vc[0],vc[1],vc[2],ii,t))
# else:
# pass
# else:
# if
O.bodies[ii].state.pos[2]<(z1+rmax)-(O.bodies[ii].state.pos[0])*tan((pi*(90-angle))/180):
# lines=f.readlines()
# for x in lines:
#
old=(re.split('\t|\n',x))
#
yref0=float(old[1])
#
vref0=float(old[4])
# if
yref0-mi<=round(O.bodies[ii].state.pos[1],6)<=yref0+mi:
# pass
# else:
# if
round(O.bodies[ii].state.vel[1],2)==vref0:
# pass
# else:
#
cp=O.bodies[ii].state.pos
#
vc=O.bodies[ii].state.vel
# t=O.time
#
f.write("{:.6f}\t{:.6f}\t{:.6f}\t{:.2f}\t{:.2f}\t{:.2f}\t{}\t{}\n".format(cp[0],cp[1],cp[2],vc[0],vc[1],vc[2],ii,t))
# else:
# pass
#
#
#
#periodflow=deltarecordtime
#nDoflow=nrecord
#massflow=[]
#stepmass=[]
#Vmean=[]
#NewVnorm=[]
#c=0
#def flowfunction():
# global c
# if len(massflow) == 0:
# M=sum(mass)
# Q=(M/deltarecord)
# massflow.append(Q)
# stepmass.append(M)
# Vm=np.mean(Vnorm)
# Vmean.append(Vm)
# c=len(Vnorm)
# else:
# M=sum(mass)
# ii=len(stepmass)-1
# Mb=stepmass[ii]
# Q=((M-Mb)/deltarecord)
# massflow.append(Q)
# stepmass.append(M)
# NewVnorm=Vnorm[c:]
# Vm=np.mean(NewVnorm)
# Vmean.append(Vm)
# c=len(Vnorm)
#
#startflowprint=(startrecordtime+(nrecord*(deltarecordtime)))
#def flowprint():
# with
open("/home/michelangelo/Scrivania/YADE_PIONA/NOZZLE_FLOW_EDGE_37/RESULTS_1/massflow.txt","a")
as f:
# step=len(massflow)
# for ii in range(step):
# mf=massflow[ii]
# vm=Vmean[ii]
# f.write("{}\t{}\n".format(mf,vm))
#
#
#startend=int((dist)/(O.dt*v))
#zp=[]
#impacted=0
#def end():
# global impacted
# for ii in shot:
#
name="/home/michelangelo/Scrivania/YADE_PIONA/NOZZLE_FLOW_EDGE_37/RESULTS_3/contact"+str(ii)+".txt"
# with open(name,"a+") as f:
# if os.stat(name).st_size == 0 :
# zp.append(0)
# else:
# zp.append(1)
# zplength=len(zp)
# impacted=zp.count(1)
# percentage=100*(impacted)*(zplength)**-1
# if percentage>95:
# O.pause()
# else:
# pass
rotationtime=(((pi*box_angle)/180)/(wt))
iterrotation=int(rotationtime/O.dt)
def motion():
O.engines[4].angularVelocity=0
O.engines[5].angularVelocity=0
O.engines[10].angularVelocity=0
translationtime_x=(2*x1*(sin(pi*box_angle/180)))/vx
itertranslation_x=int(translationtime_x/O.dt)
stoptranslation_x=iterrotation+itertranslation_x
translationtime_z=(2*x1*(1-sin(pi*(90-box_angle)/180)))/vz
itertranslation_z=int(translationtime_z/O.dt)
stoptranslation_z=stoptranslation_x+itertranslation_z
def motion_1():
O.engines[6].dead=False
O.engines[7].dead=False
def motion_2():
O.engines[6].dead=True
O.engines[7].dead=True
def motion_3():
O.engines[8].dead=False
O.engines[9].dead=False
def motion_4():
O.engines[8].velocity=0
O.engines[9].velocity=0
y_rotation=x1-x1*(1-sin(pi*(90-angle)/180))
rotationtime_2=(((pi*angle)/180)/(wt))
iterrotation_2=int(rotationtime/O.dt)
def motion_5():
O.engines[11].dead=True
O.engines[12].dead=True
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(en=MatchMaker(matches=[(Ceramic,Ceramic,0.69),(Ceramic,Aluminum,0.43),(Ceramic,Nozzle,1)]))],
[Law2_ScGeom_MindlinPhys_Mindlin()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.13,label="newton"),
RotationEngine(angularVelocity=-wt,dead=False,ids=[t_id],rotateAroundZero=True,rotationAxis=(0,0,1),zeroPoint=[0,-x1,0]),
RotationEngine(angularVelocity=wt,dead=False,ids=[t_id1],rotateAroundZero=True,rotationAxis=(0,0,1),zeroPoint=[0,-x1,0]),
TranslationEngine(velocity=vx,dead=True,ids=[t_id2],translationAxis=(1,0,0)),
TranslationEngine(velocity=-vx,dead=True,ids=[t_id3],translationAxis=(1,0,0)),
TranslationEngine(velocity=-vz,dead=True,ids=[t_id2],translationAxis=(0,1,0)),
TranslationEngine(velocity=-vz,dead=True,ids=[t_id3],translationAxis=(0,1,0)),
RotationEngine(angularVelocity=-wt,dead=False,ids=[t_id4],rotateAroundZero=True,rotationAxis=(0,0,1),zeroPoint=[0,-x1,0]),
RotationEngine(angularVelocity=-wt,dead=False,ids=nozzle,rotateAroundZero=True,rotationAxis=(1,0,0),zeroPoint=[0,y_rotation,0]),
RotationEngine(angularVelocity=-wt,dead=False,ids=pack_rotation,rotateAroundZero=True,rotationAxis=(1,0,0),zeroPoint=[0,y_rotation,0]),
#PyRunner(command='txt()',iterPeriod=periodrecord,nDo=ndorecord),
#PyRunner(command='txt_3()',firstIterRun=start,iterPeriod=periodrecord),
#PyRunner(command='flowfunction()',iterPeriod=periodflow,firstIterRun=endrecordtime,nDo=nDoflow),
#PyRunner(command='flowprint()',iterPeriod=100,firstIterRun=startflowprint,nDo=1),
PyRunner(command='motion()',iterPeriod=1,firstIterRun=iterrotation,nDo=2),
PyRunner(command='motion_1()',iterPeriod=1,firstIterRun=iterrotation,nDo=2),
PyRunner(command='motion_2()',iterPeriod=1,firstIterRun=stoptranslation_x,nDo=2),
PyRunner(command='motion_3()',iterPeriod=1,firstIterRun=stoptranslation_x,nDo=2),
PyRunner(command='motion_4()',iterPeriod=1,firstIterRun=stoptranslation_z,nDo=2),
PyRunner(command='motion_5()',iterPeriod=1,firstIterRun=iterrotation_2,nDo=2),
]
O.saveTmp()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
