对非球颗粒进行建模的方法，主要分两大类：一类是，基于球形离散单元，利用球体捆绑技术逼近目标形状颗粒，即球簇离散单元模型；另一类是，基于精确几何描述方法的颗粒模型，比如多面体颗粒、超椭球颗粒等。关于非球颗粒的建模方法，可参考综述文章 [3] 。
3.1.1. 运动方程
与球形颗粒相比，非球颗粒运动的牛顿 – 欧拉方程则采用更为一般的形式，如下：
(3.1)\[m\dot{v_i} = F_i\]
式中， \(i\in \{1,2,3\}\) 表示全局笛卡尔坐标轴； \(m\) 为颗粒质量； \(v_i\) 为颗粒平动速度，顶部点号表示对时间求导，下同； \(F_i\) 为作用在颗粒质心的合外力。
(3.2)\[I_{i}\dot{\omega}_i - (I_{j} - I_{k})\omega_{j} \omega_{k}=M_{i}\]
式中， \(i \in \{1,2,3\}\) 表示惯性主轴， \(i\) 、 \(j\) 、 \(k\) 是轮换指标； \(I_i\) 为主惯性矩； \(\omega_i\) 为角速度； \(M_i\) 为绕颗粒质心的合力矩。作用在颗粒上的力包括体积力和接触力。与球形颗粒间的接触力相比，非球颗粒间接触力的一显著特征是：接触力的法向分量（法向接触力）通常不通过颗粒的质心，从而可以产生绕质心的转矩。该转矩可以促进或抑制非球颗粒的转动。与球形颗粒简化欧拉方程比较，非球颗粒的转动方程式 (3.2) 涉及到了主惯性矩的不对称性。描述颗粒运动的牛顿和欧拉方程可分别通过标准和扩展跳跃算法来求解。此外，为方便数值处理，局部笛卡尔坐标轴沿颗粒转动惯量在初始时刻的主轴方向。
此外，相比球形颗粒离散单元模型，非球颗粒离散单元模型中的接触检测要复杂得多，这也是两者之间的一个重要不同之处。已有研究表明，在使用离散单元法进行数值模拟过程中，占绝大部分的计算资源都消耗在了颗粒间的接触检测的计算上。对于非球颗粒离散单元模拟，则更是如此，并且较球形颗粒离散单元模拟更为耗时。尽管如此，鉴于颗粒形状的重要性，开展非球颗粒离散元研究具有重要意义。
3.1.2. 颗粒位置的数值求解
颗粒运动可通过以时间增量 \(\Delta t\) 为步长的类似中心差分方法求解，即跳蛙（leapfrog）策略，或Verlet策略 [2] 。
已知颗粒在当前时刻 \(t\) 时的位置矢量 \(\mathbf{x}^{t}\) ，线加速度 \(\mathbf{a}^t\) 可由牛顿方程式 (3.1) 计算得出，即：
(3.3)\[\mathbf{a}^t = \frac{\mathbf{F}}{m}\]
式中， \(\mathbf{F}\) 为包含阻尼力的合外力。
分别对上一时步和下一时步的位置矢量 \(\mathbf{x}^{t-\Delta t}\) 与 \(\mathbf{x}^{t+\Delta t}\) 进行泰勒级数展开：
\[\mathbf{x}^{t-\Delta t} = \mathbf{x}^{t}-\Delta t \big (\frac{d\mathbf{x}}{dt} \big )^t + \frac{\Delta t^2}{2!} \big (\frac{d^2\mathbf{x}}{dt^2} \big )^t -\frac{\Delta t^3}{3!} \big (\frac{d^3\mathbf{x}}{dt^3} \big )^t +\dots+ (-1)^n\frac{\Delta t^n}{n!} \big (\frac{d^n\mathbf{x}}{dt^n} \big )^t + O(\Delta t^{n+1})\]
\[\mathbf{x}^{t+\Delta t} = \mathbf{x}^{t}+\Delta t \big (\frac{d\mathbf{x}}{dt} \big )^t + \frac{\Delta t^2}{2!} \big (\frac{d^2\mathbf{x}}{dt^2} \big )^t +\frac{\Delta t^3}{3!} \big (\frac{d^3\mathbf{x}}{dt^3} \big )^t +\dots+ \frac{\Delta t^n}{n!} \big (\frac{d^n\mathbf{x}}{dt^n} \big )^t + O(\Delta t^{n+1})\]
因此， \(\mathbf{a}^t\) 也可写成位置矢量关于时间步中心差分的形式，即：
(3.4)\[\mathbf{a}^t = \frac{\mathbf{x}^{t-\Delta t}-2\mathbf{x}^{t}+\mathbf{x}^{t+\Delta t}}{\Delta t^2} + O(\Delta t^2)\]
舍去截断误差，整理后得到：
(3.5)\[\mathbf{x}^{t+\Delta t} = 2\mathbf{x}^{t}-\mathbf{x}^{t-\Delta t}+\mathbf{a}^t\Delta t^2=\mathbf{x}^{t}+\Delta t\big (\frac{\mathbf{x}^{t}-\mathbf{x}^{t-\Delta t}}{\Delta t} + \mathbf{a}^t\Delta t \big )\]
在时刻 \(t-\frac{\Delta t}{2}\) 时的平动速度 \(\mathbf{v}^{t-\frac{\Delta t}{2}}\) 可由当前时步与前一时步的位置矢量求得，即：
(3.6)\[\mathbf{v}^{t-\frac{\Delta t}{2}} = \frac{\mathbf{x}^{t}-\mathbf{x}^{t-\Delta t}}{\Delta t}\]
在时刻 \(t+\frac{\Delta t}{2}\) 时的平动速度 \(\mathbf{v}^{t+\frac{\Delta t}{2}}\) 可由下式得出：
(3.7)\[\mathbf{v}^{t+\frac{\Delta t}{2}} = \mathbf{v}^{t-\frac{\Delta t}{2}}+\mathbf{a}^t \Delta t\]
将式 (3.6) 代入式 (3.5) 得：
(3.8)\[\mathbf{x}^{t+\Delta t} = \mathbf{x}^{t}+\Delta t\big (\mathbf{v}^{t-\frac{\Delta t}{2}} + \mathbf{a}^t\Delta t \big )\]
将式 (3.7) 代入式 (3.8) 得：
(3.9)\[\mathbf{x}^{t+\Delta t} = \mathbf{x}^{t}+\mathbf{v}^{t+\frac{\Delta t}{2}}\Delta t\]
综上，由式 (3.3) 、 (3.7) 以及 (3.9) 便可计算每一时步的位置矢量。
3.1.3. 颗粒方向的数值求解
在颗粒旋转的计算中，转动可用绕某一特定轴（欧拉轴，以单位向量 \(\mathbf{e}\) 表示）旋转一定角度 \(\theta\) 表示，即轴
— 角表示法。因此，每一时步中的颗粒旋转向量可表示为：
\[\mathbf{\theta} = \theta \mathbf{e}\]
对于球形颗粒，主惯性矩满足 \(I_{1}=I_{2}=I_{3}=I\) ，当前时步的角加速度 \(\mathbf{\beta}^{t}\) 由欧拉方程式 (3.2) 计算得到，即：
\[\mathbf{\beta} = \frac{\mathbf{M}}{I}\]
式中， \(\mathbf{M}\) 为考虑了阻尼力后的合力矩。
类似式 (3.7) ，我们可以得到在时刻 \(t+\frac{\Delta t}{2}\) 时的角速度 \(\mathbf{\omega}^{t+\frac{\Delta t}{2}}\) ，即：
(3.10)\[\mathbf{\omega}^{t+\frac{\Delta t}{2}} = \mathbf{\omega}^{t-\frac{\Delta t}{2}}+\mathbf{\beta}^t \Delta t\]
因此，旋转向量 \(\mathbf{\theta}^{t+\Delta t}\) 为
\[\mathbf{\theta}^{t+\Delta t} = \mathbf{\omega}^{t+\frac{\Delta t}{2}}\Delta t\]
对比欧拉角表示法，利用四元数 \(q\) 表示颗粒旋转能带来更多数值求解上的方便。
\[q(q_w,q_x,q_y,q_z) = e^{\frac{\theta}{2}(e_x \mathbf{i}+e_y \mathbf{j}+e_z\mathbf{k})} = \cos\frac{\theta}{2}+(e_x \mathbf{i}+e_y \mathbf{j}+e_z\mathbf{k})\sin\frac{\theta}{2}\]
下一时步的旋转角 \(\theta^{t+\Delta t}\) 和旋转轴 \(\mathbf{e}^{t+\Delta t}\) 分别为：
\[\theta^{t+\Delta t} = |\mathbf{\theta}^{t+\Delta t}|\]
\[\mathbf{e}^{t+\Delta t} = \frac{\mathbf{\theta}^{t+\Delta t}}{\theta^{t+\Delta t}}\]
四元数增量 \(\Delta q^{t+\Delta t}\) 为：
\[\Delta q^{t+\Delta t} = \cos\frac{\theta^{t+\Delta t}}{2}+(e_x^{t+\Delta t} \mathbf{i}+e_y^{t+\Delta t} \mathbf{j}+e_z^{t+\Delta t}\mathbf{k})\sin\frac{\theta^{t+\Delta t}}{2}\]
下一时步的颗粒方向（或总旋转量）以四元数表示为：
\[q^{t+\Delta t} = \Delta q^{t+\Delta t} q^{t}\]
对于非球颗粒，由于惯性矩张量 \(\mathbf{I}\) （见式 (3.11) ）中的三个主惯性矩不相等，只能通过求解一般形式的欧拉方程（见式 (3.2) ）得到颗粒方向。
(3.11)\[\begin{split}\mathbf{I}=\begin{bmatrix}
I_1 & 0 & 0 \\
0 & I_2 & 0 \\
0 & 0 & I_3
\end{bmatrix}\end{split}\]
因此，颗粒方向的求解不能像球形颗粒那样通过简单的蛙跳算法计算。此处采用一种所谓的扩展蛙跳算法 [1] ，数值求解过程简述如下：
给定在时刻 \(t-\frac{\Delta t}{2}\) 时的角动量 \(L^{t-\frac{\Delta t}{2}}\) 和时刻 \(t\) 时的外力矩 \(M^{t}\) ，在 \(t\) 时刻以及 \(t+\frac{\Delta t}{2}\) 时刻的角动量为：
\[L^{t} = L^{t-\frac{\Delta t}{2}} +M^{t}\frac{\Delta t}{2}\]
\[L^{t+\frac{\Delta t}{2}} = L^{t-\frac{\Delta t}{2}} +M^{t}\Delta t\]
因此，对应的颗粒局部坐标系下的角速度为：
\[\hat{\omega}^{t} = \mathbf{I}^{-1}TL^{t}\]
\[\hat{\omega}^{t+\frac{\Delta t}{2}} = \mathbf{I}^{-1}TL^{t+\frac{\Delta t}{2}}\]
式中， \(T\) 为在时刻 \(t\) 时由全局坐标向局部坐标变换的旋转矩阵，可通过对应的 \(q^{t}\) 求得。
在时刻 \(t\) 时， \(q^{t}\) 对时间的变化率 \(\dot{q}^{t}\) 可由下式得出：
\[\begin{split}\begin{bmatrix}
\dot{q}^{t}_w \\
\dot{q}^{t}_x \\
\dot{q}^{t}_y \\
\dot{q}^{t}_z
\end{bmatrix}=\frac{1}{2}\begin{bmatrix}
q^{t}_w & -q^{t}_x & -q^{t}_y & -q^{t}_z \\
q^{t}_x & q^{t}_w & -q^{t}_z & q^{t}_y \\
q^{t}_y & q^{t}_z & q^{t}_w & -q^{t}_x \\
q^{t}_z & -q^{t}_y & q^{t}_x & q^{t}_w
\end{bmatrix} \begin{bmatrix}
0 \\
\hat{\omega}^{t}_x \\
\hat{\omega}^{t}_y \\
\hat{\omega}^{t}_z
\end{bmatrix}\end{split}\]
\[q^{t+\frac{\Delta t}{2}} = q^{t} +\dot{q}^{t}\frac{\Delta t}{2}\]
同理可得到在时刻 \(t+\frac{\Delta t}{2}\) 时的 \(\dot{q}^{t+\frac{\Delta t}{2}}\) 。
因此，可以得到下一时步的颗粒方向 \(q^{t+\Delta t}\) ：
\[q^{t+\Delta t} = q^{t} + \dot{q}^{t+\frac{\Delta t}{2}} \frac{\Delta t}{2}\]
以及在时刻 \(t+\frac{\Delta t}{2}\) 时的全局角速度 \(\omega^{t+\frac{\Delta t}{2}}\) ：
\[\omega^{t+\frac{\Delta t}{2}} = T^{-1} \hat{\omega}^{t+\frac{\Delta t}{2}}\]