Simple model which has balls in box.
Balls are pressed by the wall at the top of the box using servo mechanism.
The target force exerted on the wall is 100N.
Wall contact force, velocity and displacement are recorded and plotted by Python and PFC3D.
The analysis is performed in deterministic mode.
Question 1 :
However, the two methods gives quite different histories for wall velocity and contact force.
This is plot from Python
This is plot from PFC3D
Why the two methods give different results and what is correct one?
Question 2 :
According to the sign convention, wall contact force must be positive for this test. However, history of wall contact force obtained from PFC3D command gives negative value. Why?
The wall velocity history is quite different. Wall velocity history from PFC3D history shows +1.0m/s. It means that the wall is moving upward. This is very confusing results.
Question 3 :
I think that the wall velocity should approach to zero when it reached the target contact force. However, the measured wall velocity still remains maximum velocity limit even when the displacement of the wall reached stationary condition. Is that the wall velocity history is not real measurement?
Results for increased target force (10000N) are as follows.
Contact force (10500N vs 8200N) and velocity histories are still different.
"A.py"
get_ipython().magic('reset -sf') #command to initialize iPython environment
get_ipython().magic('clear') #iPython screen clear
import itasca as it
import numpy as np
import config as cfg
import importlib
importlib.reload(cfg)
it.command("python-reset-state false") #the model new command will not affect python environment
it.command("""
model new
model configure dynamic
model deterministic off
model large-strain on
model title 'test code'
""")
it.set_threads(96)
print('Program start')
print('No. of threads = ', it.threads())
# SET MODEL DOMAIN
it.set_domain_min(cfg.g_domain_min)
it.set_domain_max(cfg.g_domain_max)
# BALL PROPERTIES
cfg.g_C_MODEL = 'linear' # linear hertz ...
cfg.g_ball_prop_linear = {'density':2700, 'kn':1.0e6, 'ks':1.0e3, 'fric':0.2, 'dp_n':0.7, 'dp_s':0.1}
cfg.g_facet_prop_linear = {'kn':1.0e6, 'ks':1.0e3, 'fric':0.2, 'dp_n':0.7, 'dp_s':0.1}
it.command(f"""
contact cmat default type ball-ball model {cfg.g_C_MODEL} property kn {cfg.g_ball_prop_linear['kn']} \
ks {cfg.g_ball_prop_linear['ks']} fric {cfg.g_ball_prop_linear['fric']} dp_mode 3 dp_nratio {cfg.g_ball_prop_linear['dp_n']} dp_sratio {cfg.g_ball_prop_linear['dp_s']}
contact cmat default type ball-facet model {cfg.g_C_MODEL} property kn {cfg.g_facet_prop_linear['kn']} \
ks {cfg.g_facet_prop_linear['ks']} fric {cfg.g_facet_prop_linear['fric']} dp_nratio {cfg.g_facet_prop_linear['dp_n']} dp_sratio {cfg.g_facet_prop_linear['dp_s']}
""")
# CREATE WALLS
def wall_gen(n, name, p1, p2, p3, p4):
it.command(f"""
wall generate id {n} name '{name}' polygon {p1} {p2} {p3} {p4}
""")
p1, p2, p3, p4 = (-0.1, -0.1, -0.05), (0.1, -0.1, -0.05), (0.1, 0.1, -0.05), (-0.1, 0.1, -0.05)
wall_gen(1, 'bottom', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.1, -0.1, 0.05), (0.1, -0.1, 0.05), (0.1, 0.1, 0.05), (-0.1, 0.1, 0.05)
wall_gen(2, 'top', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, -0.1, -0.05), (-0.08, 0.1, -0.05), (-0.08, 0.1, 0.05), (-0.08, -0.1, 0.05)
wall_gen(3, 'side-1', p1, p2, p3, p4)
p1, p2, p3, p4 = (0.08, -0.1, -0.05), (0.08, 0.1, -0.05), (0.08, 0.1, 0.05), (0.08, -0.1, 0.05)
wall_gen(4, 'side-2', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, -0.08, -0.05), (0.08, -0.08, -0.05), (-0.08, -0.08, 0.05), (0.08, -0.08, 0.05)
wall_gen(5, 'front', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, 0.08, -0.05), (0.08, 0.08, -0.05), (-0.08, 0.08, 0.05), (0.08, 0.08, 0.05)
wall_gen(6, 'back', p1, p2, p3, p4)
# CREATE BALLS
def create_ball(xy_lim, z_lim, r_min, r_max, n_balls):
it.command(f"""
ball generate box {-xy_lim} {xy_lim} {-xy_lim} {xy_lim} {-z_lim} {z_lim} number {n_balls} radius {r_min} {r_max}
""")
create_ball(0.08, 0.05, 0.001, 0.005, 10000)
it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model gravity 9.81
model solve ratio 1e-3
""")
create_ball(0.08, 0.05, 0.001, 0.005, 10000)
it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")
create_ball(0.08, 0.05, 0.001, 0.005, 10000)
it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")
create_ball(0.08, 0.05, 0.001, 0.005, 10000)
it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")
it.add_save_variable('cfg')
it.command(f"""
model save 'gen_balls.sav'
""")
"B.py"
get_ipython().magic('reset -sf') #command to initialize iPython environment
import itasca as it
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
it.command("python-reset-state false") #the model new command will not affect python environment
# MODEL CONFIGURATION
it.command("""model new
""")
it.command(f"""
model restore 'gen_balls.sav'
""")
target_force = (0, 0, -100)
top_cap = it.wall.find(2)
top_cap.servo_set_force(target_force)
top_cap.servo_set_vmax(0.1)
top_cap.servo_set_active(True)
print('Target Force (N) is ', top_cap.servo_force())
# HISTORY PARAMETERS
g_ex_1, g_ex_2, g_ex_3, g_ex_4, g_ex_5 = [], [], [], [], []
g_ex_2 = []
g_ex_3 = []
g_ex_4 = []
g_ex_5 = []
#g_ex_6 = []
#g_ex_7 = []
plt.close()
at = np.array([])
at = it.ballarray.pos()
print('Max. z position', np.max(at[:,2]))
max_ball_id = np.argmax(at[:,2]) + 1
print('Max. ball id', max_ball_id) # id
max_ball = it.ball.find(max_ball_id)
print('Max ball position in vector', max_ball.pos())
radius = max_ball.radius()
print('Radius', radius)
target_position = np.max(at[:,2]) + radius * 3
print('Wall target position', target_position)
top_cap.set_pos_z(target_position)
fig = plt.figure(figsize=(15, 6))
gs = gridspec.GridSpec(nrows=2, ncols=2, height_ratios=[1, 1], width_ratios=[1, 1])
def draw_plot(*args):
global g_ex_1, g_ex_2, g_ex_3
if it.cycle() % 1000 : return
g_ex_1.append(it.mech_age())
g_ex_2.append(top_cap.force_contact_z())
g_ex_3.append(top_cap.disp_z())
g_ex_4.append(top_cap.vel_z())
g_ex_5.append(top_cap.servo_gain())
plt.figure(1)
sub_plot_1 = plt.subplot(gs[0])
plt.xlabel('Time (s)')
plt.ylabel('Contact Force (N)')
plt.plot(g_ex_1, g_ex_2, linestyle = "-", color = "Red", linewidth = "0.5")
plt.grid(True, color = 'Gray', linestyle = '--')
sub_plot_2 = plt.subplot(gs[1])
plt.xlabel('Time (s)')
plt.ylabel('Top cap disp. (m)')
plt.plot(g_ex_1, g_ex_3, linestyle = "-", color = "Black", linewidth = "0.5")
sub_plot_2 = plt.subplot(gs[2])
plt.xlabel('Time (s)')
plt.ylabel('Top cap vel. (m/s)')
plt.plot(g_ex_1, g_ex_4, linestyle = "-", color = "Black", linewidth = "0.5")
sub_plot_2 = plt.subplot(gs[3])
plt.xlabel('Time (s)')
plt.ylabel('Gain')
plt.plot(g_ex_1, g_ex_5, linestyle = "-", color = "Black", linewidth = "0.5")
plt.tight_layout()
plt.draw()
it.set_callback("draw_plot",-10.9)
print('Servo gain factor is ', top_cap.servo_gainfactor())
print('Servo gain update is ', top_cap.servo_gainupdate())
it.command(f"""
model deterministic on
history interval 1000
model history name 'd-time' dynamic time-total
wall history name 'z-vel' velocity-z id 2
wall history name 'z-force' force-contact-z id 2
model solve time 0.5
""")
"Config.py"
import numpy as np
g_domain_max = (0.1, 0.1, 0.1)
g_domain_min = (-0.1, -0.1, -0.1)
# BALL PROPERTIES
g_C_MODEL = '' # linear hertz ...
g_ball_prop_linear = {}
g_facet_prop_linear = {}
