Hi all
I'm a new user of LIGGGHTS 2.2. I'm trying to build a simulation to calculate slury flow with large particles.
Unfortunately I have one problem. I need to change the drag force on the particle with the velocity. (Ex. drag = 0.0 when v = 0.0, drag = 1.0 when v > 0.0 ) So I try to use the if-statement in the script but I can not access to the velocity information of the particle in the if statement. I don' why. Maybe I guess if-statement part is no good.
Please give me advise and helpful information !
# Slary Flow in Inclined Pipe
atom_style granular
atom_modify map array
boundary m m p
newton off
echo both
communicate single vel yes
units si
region reg block -0.05 0.05 -0.05 0.05 0. 0.15 units box
create_box 1 reg
neighbor 0.02 bin
neigh_modify delay 0
# Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 1e8
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
# New pair style
pair_style gran model hertz tangential history
pair_coeff * *
variable TimeStep equal 0.00001
timestep ${TimeStep}
# New variable
variable C11 equal 1.0
variable rho equal 2500.0
variable PI equal 3.14159
variable G equal 9.81
variable mu equal 1.0
variable dP equal 1e3
variable A equal 0.05
variable alpha equal 0
variable gx equal 20.0
variable gy equal -sin(${alpha}*${PI}/180)
variable gz equal -cos(${alpha}*${PI}/180)
# External flow
variable vel1 atom 1.0/(4.0*${mu})*(${dP})*(${A}^2-(x^2+y^2))
# Froce on paraticles
variable vel0 atom sqrt(vx^2+vy^2+vz^2)
variable dragX atom -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0-v_vel1)*(v_vel0-v_vel1)*vx/v_vel0
variable dragY atom -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0-v_vel1)*(v_vel0-v_vel1)*vy/v_vel0
variable dragZ atom -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0-v_vel1)*(v_vel0-v_vel1)*vz/v_vel0
variable buoyX atom -${rho}*${G}*(4.0/3.0*${PI}*r^3)*${gx}
variable buoyY atom -${rho}*${G}*(4.0/3.0*${PI}*r^3)*${gy}
variable buoyZ atom -${rho}*${G}*(4.0/3.0*${PI}*r^3)*${gz}
label loop
variable i loop 500
variable index equal v_vel1[$i]
if "${index} <= 1e-5" then &
"variable dragX[$i] equal 0.0" &
"variable dragY[$i] equal 0.0" &
"variable dragZ[$i] equal ${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel1[$i])*(v_vel1[$i])" &
else &
"variable dragX[$i] equal -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0[$i]-v_vel1[$i])*(v_vel0[$i]-v_vel1[$i])*vx[$i]/v_vel0[$i]" &
"variable dragY[$i] equal -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0[$i]-v_vel1[$i])*(v_vel0[$i]-v_vel1[$i])*vy[$i]/v_vel0[$i]" &
"variable dragZ[$i] equal -${C11}*${rho}*(2.0*r^2*${PI})*abs(v_vel0[$i]-v_vel1[$i])*(v_vel0[$i]-v_vel1[$i])*vz[$i]/v_vel0[$i]"
next i
jump in.test loop
fix gravi all gravity ${G} vector ${gx} ${gy} ${gz}
fix flud0 all addforce v_dragX v_dragY v_dragZ
fix flud1 all addforce v_buoyX v_buoyY v_buoyZ
# Box walls
fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.05 0. 0.
# Ddstributions for insertion
fix pts1 all particletemplate/sphere 5330 atom_type 1 density constant 2500 radius constant 0.002
fix pts2 all particletemplate/sphere 5330 atom_type 1 density constant 2500 radius constant 0.003
fix pdd1 all particledistribution/discrete 1 2 pts1 0.7 pts2 0.3
# Region and insertion
group nve_group region reg
region bc cylinder z 0.0 0.0 0.045 0.00 0.0603 units box
fix ins nve_group insert/pack seed 1 distributiontemplate pdd1 &
maxattempt 100 insert_every 200000 overlapcheck yes all_in yes vel constant 0.0 0.0 0. &
region bc particles_in_region 10 ntry_mc 10000
# Apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere
# Insert the first particles so that dump is not empty
run 1
dump dmp all custom 500 post/dump*.slary id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
# Insert particles
run 200000 upto
unfix ins
Best & Thanks.
Kazuhisa