4.2.1.1废石碰撞运动分析

废石碰撞过程是种复杂且不确定的运动，其过程分析极为困难，为简化分析，本文假设废石充填颗粒为刚性球体，废石间发生的碰撞均为弹性碰撞，碰撞无粘附和滑移，利用碰撞恢复系数反应废石运动过程中由于碰撞导致的动能损失，且不考虑由碰撞产生的旋转等其他形式。

以第一次完整碰撞为例（图4-5），废石首先以入射速度与坡面碰撞V0，考虑碰撞能量损失后以出射速度V0-1反弹做抛物线运动后以经过重力加速度作用后的新的入射速度V1进行下一次碰撞。

图4-5 碰撞过程示意图

充填料堆面角α，废石颗粒碰撞入射速度V0，其水平和竖直分量为V0h和V0v，相对于坡面的切向和法向分量为V0n和V0t；碰撞后出射速度水平和竖直分量分别为V0-1h和V0-1v，切向和法向分量分别是V0-1n和V0-1t（V0v与V0h见4.2.1.2节）。首次入射速度相对于坡面的切向和法向分量V0n和V0t为公式（4-1）：

（4-1）

假设同一废石每次碰撞恢复系数相同，并用法向恢复系数e n和切向恢复系数e t来表征碰撞后的速度恢复程度，碰撞后的出射速度为公式（4-2）：

（4-2）

废石颗粒碰撞后沿堆面抛物线运动，在坡面法向方向上，根据抛物线运动的独立性原则，其运动受中立加速度在法向上的分量（gcosα）影响，法向上的速度具有对称性，即碰撞后的出射速度的法向分量与下一次碰撞后的入射速度大小相等，仅方向相反。同时根据法向上的运动情况可以求得废石抛物线运动的时间为公式（4-3）：

（4-3）

而对于沿坡面切向法向的速度分量，其速度会在沿坡面切向方向的重力加速度（gsinα）加速t秒，因此可得速度在切向上的分量为公式（4-4）：

（4-4）

于是，根据上述分析结合公式（4-3）与公式（4-4）可以得到第二次碰撞时的入射速度V1相对于坡面的切向分量V1t和法向分量V1n为公式（4-5）：

（4-5）

将公式（4-2）代入公式（4-5）可得公式（4-6）：

（4-6）

将公式（4-1）代入公式（4-6）可以得到第一次碰撞后废石切向速度分量V1t与法向速度分量V2n，以此类推，m次碰撞后废石颗粒切向速度分量Vmt和法向速度分量Vmn分别为公式（4-7）：

（4-7）

结合废石碰撞过程的力学分析[36]与碰撞恢复系数计算模型[72]，并在500 mm废石对应系数表（表4-1）的基础上可以得到不同粒径的废石的碰撞恢复系数为公式（4-8）：

表4-1 500 mm废石粒径的法向恢复系数和切向恢复系数表

坡面特征		较完整岩体	密实堆积岩块	采场壁面	松散堆积岩体
法向恢复系数en		0.7	0.5	0.5	0.3
切向恢复系数et		0.87~0.92	0.83~0.87	0.8~0.83	0.78~0.82
						（4-8）

式中：d-废石粒径，m。

4.2.1.2废石充填井抛投过程废石运动分析

结合充填井出口处废石流动特征与废石间相互碰撞的能量损耗[28]，可得出口处废石速度V为公式（4-9）：

（4-9）

式中：e-能量剩余系数，0.96；

   h1-充填井垂直高度，m；

   dmax-废石最大粒径，m。

因此，首次首次入射颗粒垂直分量V0v与水平分量V0h为公式（4-10）：

（4-10）

出废石充填井后，废石在采场内做抛物线运动。若在此运动中未与采场壁碰撞，则废石颗粒落地点到溜井出口处的水平距离S0为公式（4-11）：

（4-11）

式中：h2-采场高度，m；

4.2.1.3坡面运动中废石运动分析

废石落入堆积体后的运动过程包括碰撞和滚动两部分。废石首先沿着堆积体坡面向下碰撞弹跳，当废石垂直坡面速度降为零后，废石不再碰撞弹跳，转而发展为沿坡面向下滚动，当废石颗粒动能受坡面摩擦阻力减小至0时，废石颗粒不再滚动。

当废石颗粒垂直坡面分速度降为零，颗粒将沿坡面滚动。废石颗粒运行的距离由碰撞之间位移和碰撞结束后滚动位移组成。废石颗粒第m次碰撞与第m+1次碰撞过程中运动的距离为公式（4-12）：

（4-12）

当碰撞次数为n+1时，碰撞过程中废石颗粒运行的总距离为公式（4-13）：

（4-13）

废石颗粒滚动距离可得公式（4-14）：

（4-14）

式中：μ-废石滚动外摩擦系数，0.58。

则废石颗粒从废石充填井口到最终停止位置运行的总水平距离可得为公式（4-15）：.

（4-15）

4.2.2.1模型参数

假设充填井角度α为90、充填管高度h1为6 m、充填采场高度为h2为60 m。废石碰撞恢复系数主要受坡面体力学性质和废石尺寸影响，根据相关研究可知废石在第3次碰撞后相对于坡面的切向速度分量仅为第一次碰撞的2.7%，因此本文将碰撞次数设置为3。

根据3.1.2.1节粒级组成测试可得废石粒级组成为表4-2。因此，在模拟中可以利用1~100之间的随机数确定模型中废石粒径参数d的大小，即先采用各粒级所占百分比作为概率确定废石粒径d所处的粒级范围，然后在这一粒级范围内等概率随机生成最终的废石粒径d。例如，当生成的随机数为66.66时，此时根据表4-2可知废石对应粒径范围为200~300 mm，于是再根据[200,300]这一区间以等概率方式随机生成小数，从而最终确定废石粒径d。

假设采场尺寸长×宽×高=20 m×12 m×20 m，该立方体体积为20×12×20 m=4800 m3。根据3.1.2.3节容重测试得到废石的松散容重1616.25 kg/m3，可计算出将该立方体最大可容纳废石约7758 t。

为确定模拟中需测试的废石数，本文以各粒径均值的一半作为半径所构成的球体体积近似为废石体积（+300 mm粒径假设以350/2 mm作为该粒径范围的废石近似球体的半径），乘以3.1.2.4节密度测试所测的废石密度为2.67 t/m3，从而得到每块废石的质量。综上，根据公式（4-16）可以推出需要以表4-2对应粒级组成的废石总数n约为7449块。

（4-16）

其中：n-废石总数；

1/109为将1 mm3转换为1 m3。

表4-2 废石粒级组成

4.2.2.2模拟方案

本节以4.2.1节废石运动理论模型为基础，采用Python模拟废石偏析导致的废石粒径不均匀分布规律，并利用Matplotlib绘图库画出废石运动模型的模拟结果。具体模拟方案如下：

1. 以4.2.2.1节的方式随机生成每块废石粒径d；
2. 以废石充填井为圆心，将圆作360份等分成360个方向，假设废石运动方法依次为1°、2°、3°、......、360°、1°、2°、3°、......；
3. 对每个方向依次以随机粒径d作为参数代入公式4-15得该方向上的废石运动距离distance，以第n次模拟为例，可以得到该废石的坐标为(n%360，distance），对应粒径为d。根据4.2.2.1节可知需要废石7449块，故重复该步骤7449次，用于模拟7449块废石下放后的规律；
4. 对应每块废石可以将其在空间中的极坐标与粒径大小表示为(n%360，distance，d），以充填井作为圆心将每块废石的极坐标转换为直角坐标，可以得到直角坐标系下的废石（cos(n%360)*distance，sin(n%360)*distance，d）；
5. 对步骤4）中的7449个点采用Scipy库中的三维散点插值函数griddata进行插值，然后采用Matplotlib库中的plot_surface函数绘制出基于废石运动理论模型的模拟结果。

4.2.2.3模拟结果

（1）废石粒径总体分布规律

Matplotlib库中的plot_surface函数对废石运动理论模型的可视化模拟结果见图4-6，从该图可以看出，当充填井以与采场倾角一致的角度布置在采场中央时，废石在各方向上的运动距离具有对称性，总体上表现为大粒径废石运动距离远，小粒径废石运动距离相对较近。

图4-6 基于废石运动理论模型的模拟结果

（2）任一截面废石粒径分布规律

在对废石运动理论模型模拟结果的基础上，取任一截面废石粒径分布情况，可以得到图4-7。从该图可以看出，对于某一截面而言，废石水平距离随着废石粒径的增加而增加，且水平距离的增加值逐渐变陡，即水平距离与粒径的拟合函数为凹函数，拟合曲线斜率逐渐增大，粒径大小的增加对废石偏析的影响越来越大。因此，对于某金矿而言，有必要对废石粒径进行一定程度的控制，例如采用-150 mm以下的废石，从而能较好地控制废石偏析的程度。

图4-7 废石运动水平距离L与粒径d关系图

代码：

from scipy import interpolate
import pandas as pd
import numpy as np

from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
from pylab import *
from numpy import *
import matplotlib.pyplot as plt
import pandas as pd
import math
import numpy as np
import copy
plt.rcParams['axes.unicode_minus']=False  #用于解决不能显示负号的问题
mpl.rcParams['font.sans-serif'] = ['SimHei']

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import matplotlib.ticker as ticker

import matplotlib.colors 
import matplotlib.ticker


def v(d):
    return np.sqrt(2*0.96*g*h1)/dmax**2*d**2


def en(d):
    return 0.00028512*(d*1000)**(9/8)
#     return 0.3*d**(9/8)
#     return 0.3*d/dmax
#     return 0.3/0.4*d

def et(d):
    return 0.0047*(d*1000)**(5/6)
#     return 0.47*d**(5/6)
#     return 0.78*d/dmax+0.04
#     return 0.8/0.4*d

def v0v(d):
    return np.sqrt(v(d)**2*sintheta**2+2*g*h2)

def v0h(d):
    return v(d)*costheta

# v0v = np.sqrt(v(d)**2*sintheta**2+2*g*h2)
# v0h = v(d)*costheta

def vt(m):
    value = 0
    for i in range(0,m-1):
        value += en(d)**i*et(d)**(m-2-i)
    
#     print(et(d)**m*(v0v*sina+v0h*cosa)+2*et(d)*en(d)*tana*(v0v*cosa+v0h*sina)*value)
    return et(d)**m*(v0v(d)*sina+v0h(d)*cosa)+2*et(d)*en(d)*tana*(v0v(d)*cosa+v0h(d)*sina)*value

def vn(m):
    return en(d)**m*(v0v(d)*cosa+v0h(d)*sina)

def L1(m):
    value = 0
    for i in range(1,m):
#         print(i,(2*vn(i)*vt(i)/(g*cosa)+2*sina*vn(i)**2/(g*cosa**2)))
        value += (2*vn(i)*vt(i)/(g*cosa)+2*sina*vn(i)**2/(g*cosa**2))
#     print("value=",value)
    return value
        
def L2(m):
    return vt(m)**2/(2*g*miu*cosa)
 
def S0(d):
#     print(v(d),costheta)
#     print(d)
#     print(v(d)*costheta*(np.sqrt(v(d)**2*sintheta**2+2*g*h2)-v(d)*sintheta)/g)
    return v(d)*costheta*(np.sqrt(v(d)**2*sintheta**2+2*g*h2)-v(d)*sintheta)/g
    
def L(m,d):
#     print(m,d)
    return L1(m)+L2(m)+S0(d)

def myD():
    sizeRange = np.random.uniform(0,100)
    if sizeRange <= 2.23:
        return np.random.uniform(0,0.005)
    elif sizeRange <= 5.29:
        return np.random.uniform(0.005,0.01)
    elif sizeRange <= 12.73:
        return np.random.uniform(0.01,0.02)
    elif sizeRange <= 24.97:
        return np.random.uniform(0.02,0.04)
    elif sizeRange <= 39.5:
        return np.random.uniform(0.04,0.08)
    elif sizeRange <= 41.65:
        return np.random.uniform(0.08,0.15)
    elif sizeRange <= 63.61:
        return np.random.uniform(0.15,0.2)
    elif sizeRange <= 88.95:
        return np.random.uniform(0.2,0.3)
    elif sizeRange <= 100:
        return np.random.uniform(0.3,0.4)
    
    

a = 35
theta = 90
costheta = np.cos(np.pi*theta/180)
sintheta = np.sin(np.pi*theta/180)

cosa = np.cos(np.pi*a/180)
sina = np.sin(np.pi*a/180)
tana = np.tan(np.pi*a/180)
g = 9.8

miu = 0.58
h1 = 6  #充填管长度
h2 = 60
n = 3 #碰撞次数 
dmax = 0.4

np.random.seed(66666)

size = []
xList = []
yList = []

for thetaNew in range(0,7449,180):
    thetaNew = thetaNew%360
    
#     d = np.random.uniform(0.01,0.4)
    d = myD()
#     print(d)
    distanceTemp =  L(n,d)
    
    costhetaLocal = np.cos(np.pi*thetaNew/180)
    sinthetaLocal = np.sin(np.pi*thetaNew/180)
    
    x = distanceTemp * costhetaLocal
    y = distanceTemp * sinthetaLocal

#     print(thetaNew,d,distanceTemp,costheta,sintheta,x,y)
    xList.append(x)
    yList.append(y)
    size.append(d)

    

x = np.linspace(0,0.4,1000)

a = -0.625120779246529

b = 33.7129229042526

c = -168.828767985143

d = 1111.77680106694


y = a + b*x + c*x**2 + d*x**3

plt.figure(figsize=(14,8))

plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

plt.ylabel("废石运动水平距离L/m",fontsize=25)
plt.xlabel("粒径d/m",fontsize=25)


myColor=["grey","gold","darkviolet","turquoise","r","g","b","c",
               "k","darkorange","lightgreen","plum", "tan","khaki", "pink", "skyblue","lawngreen","salmon","grey","gold","darkviolet","turquoise","r","g","b","c",
               "k","darkorange","lightgreen","plum", "tan","khaki", "pink", "skyblue","lawngreen","salmon",
        "grey","gold","darkviolet","turquoise","r","g"]

plt.scatter(size,np.abs(xList),c=myColor,s=50)
plt.plot(x,y,label="拟合曲线 y=-0.63+33.71d-168.83d$^2$+1111.78d$^3$")

plt.legend(fontsize=20)
# plt.ylabel("粒径(mg/$m^3$)",fontsize=20)
plt.savefig("截面-反.jpg",dpi=1000)



#倒立：ax.view_init(-154,-160)


x = np.array(xList)
y = np.array(yList)
z = np.array(size)
        
X = np.linspace(np.array(x).min(), np.array(x).max(), 50)
Y = np.linspace(np.array(y).min(), np.array(y).max(), 50)
X1, Y1 = np.meshgrid(X, Y)
Z = interpolate.griddata((x, y), z, (X1, Y1), method='nearest',fill_value=0)

%matplotlib notebook


fig = plt.figure()  # 建立一个空间
ax = fig.add_subplot(111, projection='3d')  # 3D坐标
sp = ax.plot_surface(X1, Y1, Z,cmap='rainbow')
# sp= ax.contourf(X1,Y1,Z)

# ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1, 1, 1, 1]))
ax.zaxis.set_tick_params(pad=5)



#位置[左,下,右,上]
# fig.colorbar(sp,cax=fig.add_axes([0.1, 0.83, 0.8, 0.85]),orientation='horizontal')
fig.colorbar(sp,pad=0.15)
# ax.view_init(49,-55)
# ax.view_init(54,-58)
ax.view_init(27,134)
# plt.savefig("test.jpg",dpi=500)
# plt.xticks(np.arange(0,15,1))

plt.xlabel("距离下料口距离/m")
plt.ylabel("距离下料口距离/m")
ax.set_zlabel("废石粒径/m")

ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1, 1, 1, 0.85]))
# ax.set_zticks(np.array([0]))

plt.savefig("temp0329.jpg",dpi=500)