非球颗粒离散单元模型构建
非球颗粒离散单元模型
与球形颗粒相比，非球颗粒运动的牛顿 — 欧拉方程则采用更为一般的形式，如下： \begin{equation}\label{eqnewton_nonsph} E = mc^2 \end{equation} 式中，$i\in {1,2,3}$ 表示全局笛卡尔坐标轴；$m$ 为颗粒质量；$v_i$为颗粒平动速度，顶部点号表示对时间求导，下同；$F_i$为作用在颗粒质心的合外力。 \begin{equation}\label{eqeuler_nonsph} I_{i}\dot{\omega}i - (I{j} - I_{k})\omega_{j} \omega_{k}=M_{i} \end{equation} 式中，$i \in {1,2,3}$表示惯性主轴，$i$、$j$、$k$是轮换指标；$I_i$为主惯性矩；$\omega_i$为角速度；$M_i$为绕颗粒质心的合力矩。作用在颗粒上的力包括体积力和接触力。与球形颗粒间的接触力相比，非球颗粒间接触力的一显著特征是：接触力的法向分量（法向接触力）通常不通过颗粒的质心，从而可以产生绕质心的转矩。该转矩可以促进或抑制非球颗粒的转动。与球形颗粒简化欧拉方程比较，非球颗粒的转动方程式\eqref{eqeuler_nonsph}涉及到了主惯性矩的不对称性。描述颗粒运动的牛顿和欧拉方程可分别通过标准和扩展跳跃算法来求解。此外，为方便数值处理，局部笛卡尔坐标轴沿颗粒转动惯量在初始时刻的主轴方向。
此外，相比球形颗粒离散单元模型，非球颗粒离散单元模型中的接触检测要复杂得多，这也是两者之间的一个重要不同之处。已有研究表明，在使用离散单元法进行数值模拟过程中，占绝大部分的计算资源都消耗在了颗粒间的接触检测的计算上。对于非球颗粒离散单元模拟，则更是如此，并且较球形颗粒离散单元模拟更为耗时。尽管如此，鉴于颗粒形状的重要性，开展非球颗粒离散元研究具有重要意义。对非球颗粒进行建模的方法，主要分两大类：一类是，基于球形离散单元，利用球体捆绑技术逼近目标形状颗粒，即球簇离散单元模型；另一类是，基于精确几何描述方法的颗粒模型，比如多面体颗粒、超椭球颗粒等。本章内容安排如下：首先，对球簇离散单元模型作一简单介绍；然后，提出了多面体离散单元模型和超椭球离散单元模型，这两种模型在本文后续研究中分别用于棱角和光滑颗粒的离散元建模；最后，介绍了非球颗粒运动的数值求解，并自主开发了非球颗粒离散单元程序SudoDEM。
球簇离散单元模型
球簇颗粒
球簇颗粒是由一系列球（以下简称组成球）捆绑构建出的非球颗粒，主要分两种类型（见图1.1：一种是组成球间可以有重叠，也被称为clump模型；另一种是组成球之间无重叠，也被称为cluster模型。注，此处的重叠表示远大于DEM中颗粒正常接触的重叠。以上两种球簇模型的分类主要基于下述实践目的：clump模型可单纯用于表示非球颗粒，不涉及内部组成球间的接触力计算；cluster模型则涉及组成球间的接触力计算，即颗粒内部力可求解，该模型往往被用于考虑颗粒破碎的模拟。
两类球簇模型的二维示意图：(a) clump模型；(b) cluster模型
由球簇颗粒的定义可知，目标非球颗粒形状的逼近程度直接与组成球的个数相关。一般地，组成球的个数越多，球簇颗粒越能够逼近目标颗粒。此外，球簇颗粒中组成球的排布可以是任意的，因此，球簇模型具有很大的灵活性。理论上，球簇颗粒模型可以用于构建任意复杂形状的颗粒。但是，在实践过程中，组成球数量较大会显著地增加模拟计算耗时。鉴于此，有学者[@ashmawy2003]对非球颗粒球簇模型的构建进行了相关研究，尝试尽可能使用较少的组成球来构建颗粒。
接触计算
球簇颗粒的接触检测是基于组成球的接触检测，即通过检测两颗粒的组成球的接触来实现颗粒接触检测，从而大大降低了非球颗粒间的接触检测的复杂度。组成球间的颗粒接触检测（如图1.2{reference-type=“ref” reference=“figclumpcontact”}(b)）与球形离散单元模型中的接触检测（如图1.2{reference-type=“ref” reference=“figclumpcontact”}(a)）完全相同。但需要指出的是，为了加速接触检测，属于同一颗粒的组成球间可以不用进行接触检测（只针对clump模型）。球形颗粒间的距离当且仅当满足以下条件时两颗粒间存在接触： \begin{equation}\label{eqspherecontact} d_{ij} < R_i + R_j \end{equation} 式中，$d_{ij}$为两球形颗粒质心间的距离；$R_i$和$R_j$分别为两颗粒的半径。对于相互接触的两球形颗粒，颗粒间的接触重叠量（一般为接触嵌入深度）$\delta$则为: \begin{equation}\label{eqspheredelta} \delta = R_i + R_j - d_{ij} \end{equation} 其他接触几何量亦可容易地获取，如接触法向定义为两颗粒质心连线所在方向，此处从略。值得强调的是，此处以及本文后续章节均只涉及非粘性颗粒间的接触；对于粘性颗粒间的接触，式([eqspherecontact]{reference-type=“ref” reference=“eqspherecontact”})则需要考虑颗粒间可能存在一定的间隙，从而式([eqspheredelta]{reference-type=“ref” reference=“eqspheredelta”})中的$\delta$可以取负值，即模拟颗粒分离后的粘性力。
{#figclumpcontact
width=“12cm”}
需要指出的是，在每一DEM计算循环中需要进行大量颗粒间的接触检测，而每一颗粒只与其周围的少数（相对于系统中的大量颗粒）颗粒存在接触，亦即某一给定颗粒与绝大多数颗粒是无接触的。因此，为了避免对这些大量的非接触的颗粒对使用式([eqspherecontact]{reference-type=“ref” reference=“eqspherecontact”})进行精确的接触检测判断，通常需要先利用筛除算法（sweep and prune）[@baraff1992]排除这些大量的非接触颗粒对。包围盒算法是一类有效的筛除方法，广泛使用的包围盒有球包围盒、轴对齐包围盒以及有向包围盒等。其中，轴对齐包围盒（axis-aligned bounding box，AABB）[@bergen1997]算法在DEM颗粒接触检测中最为流行。顾名思义，AABB是包围在物体外的一个六面与全局坐标轴对齐的立方体盒子，且在大多数情况下，AABB可以比较紧凑地包围颗粒。如图1.2{reference-type=“ref” reference=“figclumpcontact”}(a)所示，其中的矩形框即为轴对齐包围盒。通过比较每个轴上的投影的长度，便可检测两个AABB之间的重叠情况。如果存在使AABB的投影不重叠的轴（如图1.2{reference-type=“ref” reference=“figclumpcontact”}(a)中颗粒A和颗粒B的包围盒在x轴上的投影不重叠），则两个AABB之间不存在重叠，这也意味着两个颗粒之间不重叠。如果两个AABB之间存在重叠（如图1.2{reference-type=“ref” reference=“figclumpcontact”}(a)中颗粒A和颗粒C、颗粒D间存在AABB的重叠），则两个颗粒之间可能存在重叠（如图1.2{reference-type=“ref” reference=“figclumpcontact”}(a)中颗粒A与颗粒C接触而颗粒A与颗粒D不接触）。轴对齐包围盒算法能够排除所有类似颗粒A、B这样的颗粒对。
接触模型
如前所述，球簇颗粒间的接触本质上是来自不同颗粒的组成球之间的接触。因此，传统球形离散单元的接触模型均可直接用于球簇离散单元模型。以接触刚度模型为例，我们可以灵活地选择线刚度模型、体积刚度模型以及其他非线性刚度模型（如简化的Hertz-Mindlin模型）。线刚度模型和体积刚度模型分别在本章后续小节中具体介绍，此处从略。
多面体离散单元模型
多面体颗粒
在代数几何学中，我们可以定义一个N面凸多面体为由N个包围面${p_i|i=1,\dots,N}$构成的凸包围集，即一组不等式集： \begin{equation}\label{eqpolydef} P={(x,y,z)|a_ix+b_iy+c_iz\leq d_i,\quad{}i=1,\dots,N} \end{equation} 其中$(a_i,b_i,c_i)$是第$i$个包围平面$p_i$的法向量，$d_i$是$p_i$到原点的距离。
多面体的基本几何参量（如表面积、体积、质心坐标、惯性矩等）可通过经典的计算几何算法求得。比如，在求解多面体的体积时，可先将多面体剖分成一系列的四面体，再对这些四面体的体积求和即可。其它几何参量的求解思路类似，此处从略。
接触检测
如图1.3{reference-type=“ref” reference=“figpolydem”}为两相互接触的多面体颗粒，其中夸大显示了颗粒间的重叠部分。多面体颗粒间的接触检测已经得到了很多关注，尤其在计算几何领域。给定一对可能相互接触的多面体颗粒，记为$P_A$和$P_B$，两多面体是否接触等价为两凸集是否存在交集。
多面体颗粒间的接触检测要比球形颗粒间的接触检测复杂得多。对于球形颗粒，根据其位置和半径即可快速获得两个颗粒的重叠，见式([eqspherecontact]{reference-type=“ref” reference=“eqspherecontact”})。然而，对于多面体颗粒，重叠区域（参见图1.3{reference-type=“ref” reference=“figpolydem”}）是不规则的，并且它具有多个接触点（例如在面面或面边接触处）。因此，在大多数情况下，计算两个多面体颗粒的重叠是比较耗时的。由于在每个时间步长中需要对大量的颗粒进行相互接触检测，我们可以把接触检测分成两个过程（近似接触检测和精确检测）来减少计算时间。
为了加速接触检测过程，我们首先采用轴对齐包围盒（AABB）算法进行粗略地接触检测。轴对齐包围盒算法简介参见本章上一小节。对于使用AABB算法不能排除的颗粒对，需要执行进一步的精确检测，从而判断两颗粒之间是否存在接触。在这种情况下，我们利用分离轴（或分离面）算法进行进一步接触检测。设想存在一个分离面$S={(x,y,z) | a_sx + b_sy+ c_sz = d_s }$使得来自多面体$P_\mathrm{A}$和$P_\mathrm{B}$的所有顶点分别在分离面的两侧。有三类平面可作为候选分离面[@elias2014]：
-
多面体$P_\mathrm{A}$的包围面；
-
多面体$P_\mathrm{B}$的包围面；
-
分别由多面体$P_\mathrm{A}$的一条边和多面体$P_\mathrm{B}$的一条边构成的平面。
如果能找到这样一个分离面，那么两多面体颗粒间不存在接触，否则，两多面体颗粒相互接触。 对于两相互接触的多面体颗粒，需要进一步计算接触处的相关几何量，如接触重叠量（嵌入深度或体积）、接触法向等，相关过程见下一小节。
接触计算
在每一接触处，需要计算相应的几何参量，如接触重叠量、接触中心以及接触法切向等。其中，接触重叠量可以是接触嵌入深度或接触重叠体积，而本模型采用接触重叠体积，这是因为下文具体的接触计算方法可以方便地给出接触重叠体积。值得指出的是，颗粒与墙（平面）的接触是颗粒间接触的退化形式，因此，此处重点介绍颗粒间的接触计算。
{#figshpv
width=“7cm”}
{#figpolyint
width=“9cm”}
对于球形颗粒，接触重叠体积（见图1.4{reference-type=“ref” reference=“figshpv”}）可通过简单的几何计算求得，如下： \begin{equation}\label{eqvsph} v_c=\frac{\pi}{3}[(3r_i-h_i)h_i^2+(3r_j-h_j)h_j^2] \end{equation} 式中，$r_i$(或$r_j$)和$h_i$(或$h_j$)分别是第$i$(或$j$)个颗粒的半径和球冠高度。其他几何参数也可以通过显式表达式给出，此处从略。
相比之下，对于多面体颗粒，颗粒间的重叠区域往往是不规则多面体，并没有明确的统一理论表达式来计算相应的几何参量。学者相继提出和发展了多种针对多面体颗粒的精确接触检测的算法，比如，公共平面（CP）算法[@cundall1988]、快速公共平面（FCP）算法[@nezami2004]以及内部潜在颗粒法[@boon2012]。这些方法能够计算多面体接触的嵌入深度，但具体算法的实现比较繁琐。因此，我们利用一种基于所谓的对偶理论[@muller1978; @gunther1991]的方法来求解重叠多面体。感兴趣的读者可在文献[@elias2014; @muller1978; @gunther1991]中找到关于这种方法的更多细节，在此只给出主要的操作步骤，如下：
-
全局坐标系向局部坐标系转换。局部坐标系的原点（以下简称为新原点）需要设置在接触重叠多面体内部。对于一个新的接触，由于在上一时步中并不存在相关信息，需要遍历两个接触多面体的所有边和面直到找到一个属于重叠的点为止；通过将该点向接触重叠内部移动一小距离便可得到新原点。对于一个已经存在的接触，来自上一时步的接触点可作为新原点。
-
由模型空间向对偶空间映射。模型空间指由式([eqpolydef]{reference-type=“ref” reference=“eqpolydef”})定义多面体的空间。在模型空间中，来自两多面体的每个包围面$p$映射到对偶空间中的一个点$p^{'}$。该映射关系如下： $$p \Rightarrow p^{'}=(\frac{a}{d},\frac{b}{d},\frac{c}{d})$$ 式中，$a$、$b$、$c$以及$d$是包围面$p$的平面方程参数。
-
在对偶空间中求解所有$p^{'}$的凸包。鉴于新原点在接触重叠内部，$p^{'}$总是存在的。在模型空间中，包围面离新原点越远，在对偶空间中对应的点将离新原点越近。因此，所有对偶点的凸包对应接触重叠。
-
由上一步获得的凸包按下述映射关系映射回模型空间： $$q^{'} \Rightarrow q=(\frac{a^{'}}{d^{'}},\frac{b^{'}}{d^{'}},\frac{c^{'}}{d^{'}})$$ 式中，$q^{'}$为在对偶空间中的凸包包围面；$a^{'}$、$b^{'}$、$c^{'}$以及$d^{'}$是$q^{'}$的平面方程参数；$q$是模型空间中接触重叠的一点。在模型空间中的所有接触重叠点$q$的凸包即为所需求解的接触重叠多面体。
根据所求得的接触重叠多面体，通过经典的计算几何算法便可很容易地求解接触方向和接触点。在此，我们引入接触交线来确定接触法向。接触交线是两接触颗粒表面的相交线，参见图1.3{reference-type=“ref” reference=“figpolydem”}或图1.5{reference-type=“ref” reference=“figpolyint”}。图1.5{reference-type=“ref” reference=“figpolyint”}显示了两接触多面体（见图1.3{reference-type=“ref” reference=“figpolydem”}）的相交多面体。针对接触交线的所有组成线段，利用最小二乘法拟合出一个拟合平面，并假设该拟合平面的法向为接触法向，与该法向正交的方向为切向方向。值得注意的是，此处的切向方向严格地说是切向方向所在平面（切平面），具体的某一接触切向与两颗粒在该切向平面的相对位移有关。对于球形颗粒，接触交线严格地位于拟合平面上，见图1.4{reference-type=“ref” reference=“figshpv”}，然而对于多面体颗粒则不然。此外，可假设重叠部分（多面体）的质心作为接触点，亦即接触力的作用点。值得指出的是，上述接触法向和接触点的定义是非常近似的，因为两多面体的重叠部分（见图1.3{reference-type=“ref” reference=“figpolydem”}）是不规则的，并且实际上存在多个接触点，比如在面面接触或面边接触处。此外，在多数情况下，接触点并不落在拟合平面上。
体刚度接触模型
接触刚度模型，作为DEM中常用的接触模型，定义了接触力$F$与相对位移$\delta$之间的弹性关系，见式([eqlaw]{reference-type=“ref” reference=“eqlaw”})。为方面描述以及建模，接触力通常被分为两正交分量：接触平面上的切向接触力$F_t$和接触平面法线方向的法向接触力$F_n$，如图1.3{reference-type=“ref” reference=“figpolydem”}所示。需要指出的是，接触平面被认为是由两多面体的相交曲线的最小二乘拟合所得[@elias2014]。 \begin{equation}\label{eqlaw} F=K\delta \end{equation} 式中，$K$为接触刚度，由接触颗粒的材料特性决定。如果接触刚度$K$在整个计算模拟过程中为常数，那么式([eqlaw]{reference-type=“ref” reference=“eqlaw”})定义了一个线弹性接触模型。本文中的接触刚度$K$由式([eqK]{reference-type=“ref” reference=“eqK”})给出，其中，上标A和B表示相接触的两个颗粒。 \begin{equation}\label{eqK} K=\frac{K^AK^B}{K^A+K^B} \end{equation}
通常，式([eqlaw]{reference-type=“ref” reference=“eqlaw”})中的接触重叠$\delta$为嵌入深度。此处，采用一个适用于任意颗粒形状的法向接触模型。如式([eqvlaw]{reference-type=“ref” reference=“eqvlaw”})所示，法向接触力与接触重叠体积相关。该法向接触模型已被其他学者[@chen2010; @nassauer2013; @nassauer2013c; @elias2014]使用。 \begin{equation}\label{eqvlaw} F_n=K_nv_c \end{equation} 式中，$K_n$[N/m$^3$]为接触法向刚度（体积刚度），由式([eqK]{reference-type=“ref” reference=“eqK”})给出；$v_c$[m$^3$]为接触重叠体积，其详细计算参见上一小节。
此外，采用最常用的切向接触模型，如下： \begin{equation}\label{eqFs} \Delta F_t=-K_t \Delta u \end{equation} 式中，$K_t$ [N/m]为接触切向刚度；$\Delta u$与$\Delta F_t$分别为当前时步切向接触位移增量和切向接触力增量。
库伦条件或滑动摩擦模型用于确定接触的滑动。当$F_t$大于$\tan\phi_\mu\cdot F_n$ （其中$\phi_\mu$为滑动摩擦角）时，接触发生滑动，此时$F_t$须减小并维持$\tan\phi_\mu\cdot F_n$。需要说明的是，鉴于滚动摩擦的影响相对较小[@mack2011]，本文模型不涉及滚动摩擦。
超椭球离散单元模型
超椭球颗粒
基本定义
超椭球面是一类特殊的闭合超二次曲面，能够表达丰富的几何形状。在局部笛卡尔坐标系下，超椭球的表面方程可以定义为[@barr1981]： \begin{equation}\label{eqsurface_super} (|\frac{x}{a}|^{\frac{2}{\epsilon_1}} + |\frac{y}{b}|^{\frac{2}{\epsilon_2}})^{\frac{\epsilon_1}{\epsilon_2}} + |\frac{z}{c}|^{\frac{2}{\epsilon_2}} = 1 \end{equation} 其中，$a$、$b$以及$c$分别是在x、y和z轴方向上的半长轴长度，$\epsilon_i (i=1,2)$是确定颗粒边缘尖锐程度的形状参数。 感兴趣的读者可在文献[@williams1992; @cleary1997efficient]中找到超椭球的其他类似定义。本文主要研究凸形状，其对应的$\epsilon_i$在$(0,2)$之间。如图1.6{reference-type=“ref” reference=“figsuperellipsoids”}所示，改变$\epsilon_i$可以得到多种形状。 特别地，$\epsilon_i \rightarrow 0$时对应立方体形状，而$\epsilon_i \rightarrow 2$时对应八面体形状。
{#figsuperellipsoids
width=“8cm”}
主要几何参量
以下给出了本模型中涉及到的超椭球的一系列重要的几何参量，如体积、主惯性矩等。 超椭球的体积由下式给出: \begin{equation}\label{eqVolume} V = 2AB(\frac{1}{2}\epsilon_1+1, \epsilon_1)B(\frac{1}{2}\epsilon_2, \frac{1}{2}\epsilon_2)\end{equation} 式中， $A$等于$abc\epsilon_1\epsilon_2$； $B(x, y)$为贝塔函数，定义为 $$B(x, y) = 2\int_0^{\frac{\pi}{2}}\mathrm{sin}^{2x-1}\phi \mathrm{cos}^{2y-1}\phi d\phi = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$$ 其中，$\Gamma(x)$为伽马函数。 超椭球的主惯性矩为： \begin{equation}\label{eqinertia} \begin{aligned} I_1 &= I_{xx} =\frac{1}{2} \rho A(b^2\beta_1 + 4c^2\beta_2) \\ I_2 &= I_{yy} = \frac{1}{2} \rho A(a^2\beta_1 + 4c^2\beta_2) \\ I_3 &= I_{zz} = \frac{1}{2} \rho A(a^2 + b^2)\beta_1 \end{aligned} \end{equation} 式中， $\rho$ 为材料密度, 参数$\beta_1$和$\beta_2$定义如下： \begin{equation}\left{ \begin{array}{l@{}l} \beta_1 = B(\frac{3}{2}\epsilon_2, \frac{1}{2}\epsilon_2)B(\frac{1}{2}\epsilon_1, 2\epsilon_1 + 1)\\ \beta_2 = B(\frac{1}{2}\epsilon_2, \frac{1}{2}\epsilon_2 + 1)B(\frac{3}{2}\epsilon_1,\epsilon_1+1) \end{array} \right.\end{equation}
超椭球主曲率
给定超椭球的参数方程为： \begin{equation} \mathbf{r}(\theta,\phi)= \begin{bmatrix} \mathrm{Sign}(\cos\theta)r_x|\cos\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2} \\ \mathrm{Sign}(\sin\theta)r_y|\sin\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2} \\ \mathrm{Sign}(\sin\phi)r_z|\sin\phi|^{\epsilon_2} \end{bmatrix} \end{equation} 且 $\theta \in [0,2\pi),\phi \in [-\frac{\pi}{2},\frac{\pi}{2}]$。其中， $r_x$、$r_y$和 $r_z$ 为局部坐标系下在主轴方向的半长轴；$\mathrm{Sign}(x)$表示符号函数，定义如下： \begin{equation}\label{eqsign} \mathrm{Sign}(x) = \begin{cases} -1 & \text{if} \quad x < 0, \\ 0 & \text{if} \quad x = 0, \\ 1 & \text{if} \quad x > 0. \end{cases}\end{equation}
因此，参数方程的一阶和二阶导数为： \begin{equation}\mathbf{r}_{\theta}= \begin{bmatrix} -\mathrm{Sign}(\cos\theta)r_x\epsilon_1|\cos\theta|^{\epsilon_1 -2}|\cos\phi|^{\epsilon_2}\cos\theta\sin\theta \\ \mathrm{Sign}(\sin\theta)r_y\epsilon_1|\sin\theta|^{\epsilon_1-2}|\cos\phi|^{\epsilon_2}\cos\theta\sin\theta \\ 0 \end{bmatrix} \end{equation}
\begin{equation} \mathbf{r}_{\phi}= \begin{bmatrix} -\mathrm{Sign}(\cos\theta)r_x\epsilon_2|\cos\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2-2}\cos\phi\sin\phi \\ -\mathrm{Sign}(\sin\theta)r_y\epsilon_2|\sin\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2-2}\cos\phi\sin\phi \\ \mathrm{Sign}(\sin\phi)r_z\epsilon_2|\sin\phi|^{\epsilon_2-2}\cos\phi\sin\phi \end{bmatrix} \end{equation}
\begin{equation} \begin{aligned} \mathbf{r}_{\theta\phi}=\frac{\sin2\theta\sin2\phi}{4} \begin{bmatrix} \mathrm{Sign}(\cos\theta)r_x\epsilon_1\epsilon_2|\cos\theta|^{\epsilon_1-2}|\cos\phi|^{\epsilon_2-2} \\ -\mathrm{Sign}(\sin\theta)r_y\epsilon_1\epsilon_2|\sin\theta|^{\epsilon_1-2}|\cos\phi|^{\epsilon_2-2} \\ 0 \end{bmatrix} \end{aligned} \end{equation}
\begin{equation} \mathbf{r}_{\theta\theta}= \begin{bmatrix} \mathrm{Sign}(\cos\theta)r_x\epsilon_1|\cos\theta|^{\epsilon_1-2}|\cos\phi|^{\epsilon_2}(\epsilon_1\sin^2\theta-1) \\ \mathrm{Sign}(\sin\theta)r_y\epsilon_1|\sin\theta|^{\epsilon_1-2}|\cos\phi|^{\epsilon_2}(\epsilon_1\cos^2\theta-1)\\ 0 \end{bmatrix} \end{equation}
\begin{equation} \mathbf{r}_{\phi\phi}= \begin{bmatrix} \mathrm{Sign}(\cos\theta)r_x\epsilon_2|\cos\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2-2}(\epsilon_2\sin^2\phi-1) \\ \mathrm{Sign}(\sin\theta)r_y\epsilon_2|\sin\theta|^{\epsilon_1}|\cos\phi|^{\epsilon_2-2}(\epsilon_2\sin^2\phi-1)\\ \mathrm{Sign}(\sin\phi)r_z\epsilon_2|\sin\phi|^{\epsilon_2-2}(\epsilon_2\cos^2\phi-1) \end{bmatrix} \end{equation}
超椭球表面的外法线单位向量为：
\begin{equation} \mathbf{n}=\frac{\mathbf{r}{\theta} \times \mathbf{r}{\phi}}{| \mathbf{r}{\theta}\times \mathbf{r}{\phi} |} \end{equation}
式中， ‘$\times$‘表示向量叉乘， ‘$|\quad|$‘表示向量模。
第一基本形式系数： \begin{equation} E=\mathbf{r}{\theta}\cdot \mathbf{r}{\theta} \quad{} F=\mathbf{r}{\theta}\cdot \mathbf{r}{\phi} \quad{} G=\mathbf{r}{\phi}\cdot \mathbf{r}{\phi} \end{equation}
其中，'$\cdot$’ 表示向量点乘。 第二基本形式系数： \begin{equation} L=\mathbf{r}{\theta\theta}\cdot \mathbf{n} \quad M=\mathbf{r}{\theta\phi}\cdot \mathbf{n}=0 \quad N=\mathbf{r}_{\phi\phi}\cdot \mathbf{n} \end{equation}
高斯曲率： $$K=\frac{LN}{EG-F^2}$$ 平均曲率： $$H=\frac{EN+GL}{2(EG-F^2)}$$ 主曲率： \begin{equation}\rho_{I}=H+\sqrt{H^2-K} \quad \rho_{II}=H-\sqrt{H^2-K}\end{equation}
几何变换函数
给定超椭球表面上一点的局部球坐标 $(\theta, \phi)$，对应的笛卡尔坐标 $(x, y, z)$ 可由如下函数 $\mathbb{F}_1$获取： \begin{equation}\label{eqF} \begin{aligned} x =& \mathrm{Sign}(\mathrm{cos}(\theta))a|\mathrm{cos}(\theta)|^{\epsilon_1} |\mathrm{cos}(\phi)|^{\epsilon_2}\\ y =& \mathrm{Sign}(\mathrm{sin}(\theta))b|\mathrm{sin}(\theta)|^{\epsilon_1} |\mathrm{cos}(\phi)|^{\epsilon_2}\\ z =& \mathrm{Sign}(\mathrm{sin}(\phi))c|\mathrm{sin}(\phi)|^{\epsilon_2} \end{aligned} \end{equation}
给定超椭球表面上一法向量 $(n_x, n_y, n_z)$，对应的局部球坐标 $(\theta, \phi)$ 可由如下函数 $\mathbb{F}_2$获取： \begin{equation}\label{eqf} \begin{aligned} \theta = \mathrm{atan2}(\mathrm{Sign}(n_y)|b n_y|^{\frac{1}{2 - \epsilon_1}}, \mathrm{Sign}(n_x)|an_x|^{\frac{1}{2 - \epsilon_1}}) \\ \phi = \mathrm{atan2}(\mathrm{Sign}(n_z)|cn_z|\mathrm{cos}(\theta)|^{2 - \epsilon_2}|^{\frac{1}{2 - \epsilon_2}}, |an_x|^{\frac{1}{2 - \epsilon_2}}) \end{aligned} \end{equation} 式中，$\mathrm{atan2}(x, y)$ 为$\frac{x}{y}$ 的反正切函数，取值区间为 $(-\pi, \pi]$； $\mathrm{Sign}(x)$ 为符号函数，定义参见式([eqsign]{reference-type=“ref” reference=“eqsign”})。
接触计算
在DEM的每个循环时步中都要更新颗粒间的接触嵌入深度和接触的方向，以便获得更新的接触力。给定两个可能接触的相邻颗粒A和B，把来自两颗粒的可能的接触点分别记为$\mathbf{p}^A$和$\mathbf{p}^B$，对应的潜在接触嵌入向量记为$\mathbf{d} = (\mathbf{p}^B-\mathbf{p}^A)$，参见图1.7{reference-type=“ref” reference=“figdetection”}(a)。根据公共法向概念[@johnson1985; @wellmann2008]，待求解的接触点须要使接触嵌入深度最小，并且须要同时满足以下约束条件：
(1) 在接触点$\mathbf{p}^A$ 和 $\mathbf{p}^B$的颗粒表面外法线单位向量 $\mathbf{n}^A$ 和 $\mathbf{n}^B$满足反向平行，并且与接触方向$\mathbf{c}$平行，即： \begin{equation}\label{eqnnc} \mathbf{n}^A = -\mathbf{n}^B = \mathbf{c}\end{equation}
(2) 潜在的接触嵌入向量$\mathbf{d}$与接触方向$\mathbf{c}$平行，即： \begin{equation}\label{eqdc} \mathbf{d} \times \mathbf{c} = \mathbf{0}\end{equation} 因此，找出目标接触点是一个最优化问题。如Wellmann等[@wellmann2008]所建议的，接触方向$\mathbf{c}$在局部球面坐标系中可用两个角度进行参数化描述，即： \begin{equation}\mathbf{c}(\alpha, \beta) = \mathrm{cos}\alpha\mathrm{cos}\beta \mathbf{i}+ \mathrm{sin}\alpha\mathrm{cos}\beta \mathbf{j}+ \mathrm{sin}\beta \mathbf{k}\end{equation} 其中，$\mathbf{i}$、$\mathbf{j}$、$\mathbf{k}$是全局笛卡尔坐标系的单位基矢量。 因此，考虑式([eqnnc]{reference-type=“ref” reference=“eqnnc”})，接触点可表示为： $$% \left{ \begin{array}{ll} \mathbf{p}^A = \mathbf{T}^{-1}_A\mathbb{F}_1^A(\mathbb{F}_2^A(\mathbf{T}_A\mathbf{c}(\alpha, \beta))) + \mathbf{s}^A\\ \mathbf{p}^B = \mathbf{T}^{-1}_B\mathbb{F}_1^B(\mathbb{F}_2^B(-\mathbf{T}_B\mathbf{c}(\alpha, \beta))) + \mathbf{s}^B \end{array} \right.$$ 其中，$\mathbf{T}$是从全局坐标系到颗粒局部坐标系的旋转矩阵，$\mathbf{T}^{-1}$是其逆矩阵； $\mathbf{s}$是颗粒在全局坐标系下的位置矢量；$\mathbb{F}_1$和$\mathbb{F}_2$是颗粒局部坐标系中的两个几何变换函数，参见式([eqF]{reference-type=“ref” reference=“eqF”}) 和([eqf]{reference-type=“ref” reference=“eqf”})。 因此，寻找嵌入深度$||\mathbf{d}||$ 等价为求解以下带两个参数的无约束优化问题： \begin{equation}\label{eqfinald} \underset{\alpha, \beta}{\mathrm{min}}||\mathbf{d}|| = \underset{\alpha, \beta}{\mathrm{min}}||\mathbf{p}^B - \mathbf{p}^A||\end{equation} 本文采用Nelder-Mead单纯形算法[@lagarias1998]来获得鲁棒性的求解过程。值得注意的是，当式([eqfinald]{reference-type=“ref” reference=“eqfinald”})达到全局最小值时， 式([eqdc]{reference-type=“ref” reference=“eqdc”})自动满足[@wellmann2008]。
{#figdetection
width=“10cm”}
接触检测
鉴于在每个时间步长都需要对大量的颗粒进行接触检测，可采用近似接触检测和精确接触检测的组合算法来加速接触检测过程。首先，采用两级近似接触检测方案。在第一级，轴对齐包围盒（AABB）算法（详见本章第二节）用于排除大部分不相互接触的颗粒对。鉴于超椭球的对称性，对于每个颗粒，使用固定尺寸的立方体AABB，即无需每一时步重新计算颗粒的AABB。然后，球形包围盒用于第二级的进一步筛除。对于剩下的可能存在相互接触的颗粒对，则需要进行精确接触检测，详见上一小节。需要强调的是，在优化求解过程中，可以在每次迭代时检查以下条件： \begin{equation}\mathbf{d} \do\cdot \mathbf{c} > 0\end{equation} 如果上述条件得到满足，那么所考察的颗粒对之间不存在接触[@wellmann2008]，参见图1.7{reference-type=“ref” reference=“figdetection”}(b)，因而，可以排除当前的颗粒对，并终止该颗粒对的进一步接触检测。
线刚度接触模型
DEM中假设颗粒是完全刚性但颗粒间可以发生重叠（所谓的软球模型）[@cundall1979]，参见图1.8{reference-type=“ref” reference=“figsuperedem”}。 因此，可以根据给定接触本构模型和当前的重叠来计算排斥力。目前，学者们针对不同研究问题提出了各种接触模型。然而，由于线性弹簧模型[@bagi2003]的简易性且对于非球颗粒可以显著提升计算效率（本文给出了详细讨论），其使用非常广泛。故本文的大部分数值模拟采用了这种模型，由下式给出： \begin{equation}\left{ \begin{array}{l@{}l} F_n=\delta K_n\\ F_t=\mathrm{min}{F_t^{'}-\Delta u K_t, \mu F_n} \end{array} \right.\end{equation} 其中，$F_n$ 和 $F_t$分别是法向和切向接触力；$K_n$和$K_t$是对应的法向和切向接触刚度；$\Delta u$是每个时间步长的接触切向位移增量；$\delta$是接触处颗粒间的相互嵌入（重叠）深度；$\mu$是颗粒间的摩擦系数； $F_t^{'}$是上一个时间步长的切向接触力。需要指出的是，库仑条件（或滑动摩擦模型）被施加到切向接触力的计算。当接触形成时，切向接触力$F_t^{'}$被初始化为零。接触刚度$K_$（$$代表$n$或$t$）设置为两个接触颗粒材料刚度的调和平均值，见式 （[equK]{reference-type=“ref” reference=“equK”}），其中上标A和B表示两个接触的颗粒。
\begin{equation}\label{equK} K_{}=\frac{2K_{}^AK_{}^B}{K_{}^A+K_{*}^B}\end{equation}
{#figsuperedem
width=“8cm”}
颗粒状态的数值求解
颗粒运动可通过以时间增量$\Delta t$为步长的类似中心差分方法求解，即跳蛙（leapfrog）策略，或Verlet策略[@rapaport2004]。
位置
已知颗粒在当前时刻$t$时的位置矢量$\mathbf{x}^{t}$，线加速度$\mathbf{a}^t$可由牛顿方程式([eqnewton_sph]{reference-type=“ref” reference=“eqnewton_sph”})计算得出，即： \begin{equation}\label{eqtransacc} \mathbf{a}^t = \frac{\mathbf{F}}{m}\end{equation} 式中，$\mathbf{F}$为包含阻尼力的合外力。
分别对上一时步和下一时步的位置矢量$\mathbf{x}^{t-\Delta t}$与$\mathbf{x}^{t+\Delta t}$进行泰勒级数展开： \begin{equation}\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})\end{equation} \begin{equation}\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})\end{equation} 因此，$\mathbf{a}^t$也可写成位置矢量关于时间步中心差分的形式，即： \begin{equation}\label{eqxtdeltato} \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)\end{equation} 舍去截断误差，整理后得到： \begin{equation}\label{eqxtdeltat} \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 )\end{equation} 在时刻$t-\frac{\Delta t}{2}$时的平动速度$\mathbf{v}^{t-\frac{\Delta t}{2}}$可由当前时步与前一时步的位置矢量求得，即： \begin{equation}\label{eqvt_deltat} \mathbf{v}^{t-\frac{\Delta t}{2}} = \frac{\mathbf{x}^{t}-\mathbf{x}^{t-\Delta t}}{\Delta t}\end{equation} 在时刻$t+\frac{\Delta t}{2}$时的平动速度$\mathbf{v}^{t+\frac{\Delta t}{2}}$可由下式得出： \begin{equation}\label{eqvtdeltat} \mathbf{v}^{t+\frac{\Delta t}{2}} = \mathbf{v}^{t-\frac{\Delta t}{2}}+\mathbf{a}^t \Delta t\end{equation} 将式([eqvt_deltat]{reference-type=“ref” reference=“eqvt_deltat”})代入式([eqxtdeltat]{reference-type=“ref” reference=“eqxtdeltat”})得： \begin{equation}\label{eqxtdeltat1} \mathbf{x}^{t+\Delta t} = \mathbf{x}^{t}+\Delta t\big (\mathbf{v}^{t-\frac{\Delta t}{2}} + \mathbf{a}^t\Delta t \big )\end{equation} 将式([eqvtdeltat]{reference-type=“ref” reference=“eqvtdeltat”})代入式([eqxtdeltat1]{reference-type=“ref” reference=“eqxtdeltat1”})得： \begin{equation}\label{eqftdeltat2} \mathbf{x}^{t+\Delta t} = \mathbf{x}^{t}+\mathbf{v}^{t+\frac{\Delta t}{2}}\Delta t\end{equation} 综上，由式([eqtransacc]{reference-type=“ref” reference=“eqtransacc”})、([eqvtdeltat]{reference-type=“ref” reference=“eqvtdeltat”})以及([eqftdeltat2]{reference-type=“ref” reference=“eqftdeltat2”})便可计算每一时步的位置矢量。
方向
在颗粒旋转的计算中，转动可用绕某一特定轴（欧拉轴，以单位向量$\mathbf{e}$表示）旋转一定角度$\theta$表示，即轴 — 角表示法。因此，每一时步中的颗粒旋转向量可表示为： \begin{equation}\mathbf{\theta} = \theta \mathbf{e}\end{equation} 对于球形颗粒，主惯性矩满足$I_{1}=I_{2}=I_{3}=I$，当前时步的角加速度$\mathbf{\beta}^{t}$由欧拉方程式([eqeuler_sph]{reference-type=“ref” reference=“eqeuler_sph”})计算得到，即： \begin{equation}\mathbf{\beta} = \frac{\mathbf{M}}{I}\end{equation} 式中，$\mathbf{M}$为考虑了阻尼力后的合力矩。
类似式([eqvtdeltat]{reference-type=“ref” reference=“eqvtdeltat”})，我们可以得到在时刻$t+\frac{\Delta t}{2}$时的角速度$\mathbf{\omega}^{t+\frac{\Delta t}{2}}$，即： \begin{equation}\label{eqangvtdeltat} \mathbf{\omega}^{t+\frac{\Delta t}{2}} = \mathbf{\omega}^{t-\frac{\Delta t}{2}}+\mathbf{\beta}^t \Delta t\end{equation} 因此，旋转向量$\mathbf{\theta}^{t+\Delta t}$为 \begin{equation}\mathbf{\theta}^{t+\Delta t} = \mathbf{\omega}^{t+\frac{\Delta t}{2}}\Delta t\end{equation} 对比欧拉角表示法，利用四元数$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}$分别为： \begin{equation}\theta^{t+\Delta t} = |\mathbf{\theta}^{t+\Delta t}|\end{equation} \begin{equation}\mathbf{e}^{t+\Delta t} = \frac{\mathbf{\theta}^{t+\Delta t}}{\theta^{t+\Delta t}}\end{equation} 四元数增量$\Delta q^{t+\Delta t}$为： \begin{equation}\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}\end{equation} 下一时步的颗粒方向（或总旋转量）以四元数表示为： $$q^{t+\Delta t} = \Delta q^{t+\Delta t} q^{t}$$
对于非球颗粒，由于惯性矩张量$\mathbf{I}$（见式（[nonsph_inertialM]{reference-type=“ref” reference=“nonsph_inertialM”}））中的三个主惯性矩不相等，只能通过求解一般形式的欧拉方程（见式([eqeuler_nonsph]{reference-type=“ref” reference=“eqeuler_nonsph”})）得到颗粒方向。 \begin{equation}\label{nonsph_inertialM} \mathbf{I}=\begin{bmatrix} I_1 & 0 & 0 \\ 0 & I_2 & 0 \\ 0 & 0 & I_3 \end{bmatrix}\end{equation} 因此，颗粒方向的求解不能像球形颗粒那样通过简单的蛙跳算法计算。此处采用一种所谓的扩展蛙跳算法[@allen1989]，数值求解过程简述如下：
给定在时刻$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$$ 因此，对应的颗粒局部坐标系下的角速度为： \begin{equation}\hat{\omega}^{t} = \mathbf{I}^{-1}TL^{t}\end{equation} \begin{equation}\hat{\omega}^{t+\frac{\Delta t}{2}} = \mathbf{I}^{-1}TL^{t+\frac{\Delta t}{2}}\end{equation} 式中，$T$为在时刻$t$时由全局坐标向局部坐标变换的旋转矩阵，可通过对应的$q^{t}$求得。 在时刻$t$时，$q^{t}$对时间的变化率$\dot{q}^{t}$可由下式得出： \begin{equation}\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{equation} $$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}}$： \begin{equation}\omega^{t+\frac{\Delta t}{2}} = T^{-1} \hat{\omega}^{t+\frac{\Delta t}{2}}\end{equation}
小结
本章介绍了三类非球颗粒离散单元模型，即球簇离散单元模型、多面体离散单元模型以及超椭球离散单元模型，并对相应模型进行了编程实现，开发了非球颗粒离散元核心程序SudoDEM。