# [Q]Python and PFC3D give different histories. Why?

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

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'
""")

print('Program start')

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

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 = {}
``````

Hello,
Here are my replies to your questions:
Question 1 - In your python function you are getting wall values at cycle location -10.9. Histories will be calculated at the end of the cycle. If you use `it.set_callback("draw_plot",41.0)` (41.0 is right after force-displacement calculations) you will get the same plots between your python plots and PFC plots.

Question 2 - This is correct. The wall is moving upward and thus creating a tensile (negative) normal force between the balls and the wall. Actually, the wall is fluctuating up and down - see the next comment.

Question 3 - The wall is actually fluctuating between the two limiting velocities of 0.1 and -0.1. You can see this is you set the history interval to 1. After running the model you provided and changing the history interval to 1, here are the results:

If you zoom in you can see each cycle the wall is changing velocity from -0.1 to 0.1 m/s. The same happens for the force. This tells me the servo mechanism probably needs to be changed.