SudoSim Blog

• 新发布
• SudoSim新闻
• 3: 教程
• 3.1: SudoDEM培训
• 3.2: SudoSimGUI简明教程
• 3.3: SudoDEM系列教程：从准备到安装
• 3.4: 开源非球离散元程序SudoDEM：一图预览
• 3.5: 非球颗粒离散单元模型构建

新发布

SudoSim新闻

3 - 教程

3.1 - SudoDEM培训

视频

3.2 - SudoSimGUI简明教程

SudoSimGUI简介

SudoSimGUI支持Windows、Linux和MacOS，有Server和Client两个版本。用户可在远程服务器上运行服务端命令行程序，采用本地客户端远程连接并进行实时可视化和调试等工作。

服务端程序

在命令行执行，

sudosim -n -s hostname -p port

其中'-n’为不采用gui运行服务端程序；‘hostname’为服务器地址，默认为localhost；‘port’为服务器端口，默认为9091。

客户端程序

1. 点击主程序即可打开GUI界面。
2. 先在设置中的高级选项进行服务器地址和端口指定，默认分别为localhost和9091。
3. 点击连接服务器按钮，成功后可像运行本地模拟程序一样执行程序。
4. 点击remote view按钮，程序将自动从远程服务端下载数据并可视化显示。
5. 按照本地设置的显示刷新频率，程序自动提取远程服务器数据实时更新。

注意事项

1. 点击连接服务器后，本地程序自动变为客户端角色，在shell内执行的命令全部发送到服务器端。
2. 为了在客户端模式查看本地程序，保留了lapp对象（表示local app），以该对象开头的命令不被发送到服务端执行。

演示

3.3 - SudoDEM系列教程：从准备到安装

SudoDEM基于Linux平台开发，提供 Ubuntu发行版（14.04～18.04）的二进制安装（压缩）包，解压即可使用，无需root权限，不向系统核心目录写入文件，完全“绿色”。项目后期会提供Singularity容器的二进制包，用于任何Linux发行版。

if 你是编程新手 or Linux新手：
  建议先尝试使用编译的二进制包
  请继续阅读本文
else:
  可直接下载源代码进行编译
  你可以快速跳过本文

如果你能非常容易地看懂上面代码块表达的意思，那么你可以继续往下看，否则，请直接拖到文章末尾点赞。

使用SudoDEM，你只需要掌握最基本的Linux系统知识和Python编程技巧。本文将介绍这些必备知识。首先，请准备Linux系统，推荐Ubuntu 14~18；你可以通过以下方案获得Linux系统：

• Windows上安装虚拟机（Windows重度用户，体验Linux）
• Windows 10开启Linux子系统（Windows重度用户，体验Linux）
• 安装与Windows并存的多系统（需要重启电脑切换系统）
• 安装单一系统，如深度Deepin（Windows轻度用户）
• 安装单一系统，如Ubuntu（Linux重度用户）

[Linux最基础命令]

下面以Ubuntu系统为例（不同Linux发行版（如CentOS）命令可能稍有不同），介绍几个最为基本的命令和快捷方式：

• 打开终端（terminal，即命令行窗口）：CTRL+ALT+T
• 安装软件到系统：sudo apt-get install package 或者sudo apt install package
• 显示当前路径（工作目录）：pwd
• 改变路径（工作目录）：cd directory

路径/目录directory可以是相对路径或绝对路径。

其中，绝对路径是以Linux根目录（/）开始，比如：/bin，/home/sway等。为了方便访问用户目录（如/home/sway/），系统给出了用户目录的缩写目录，表示为~；～就等价为前述/home/sway。注意， Linux系统下可以存在不同用户的用户目录，但~符号只表示当前登录用户（账号）的用户目录；并且，当前用户对其他用户的用户目录不具有访问权限。值得一提的是，Linux是天然多用户系统，即允许多个用户同时登录同一系统（但有各自的运行环境，这些配置保存在各自的用户目录下；注：“同时登录”不是访问同一桌面）。比如用户sway已经登录系统进行工作，用户SC则可以通过网络（如ssh远程登录命令）进入系统工作；不同用户互不干扰（除了抢占硬件资源，比如大家都在用SudoDEM跑计算）。 相对路径就是基于某一特定目录。严格意义上讲，基于根目录的路径也可以叫做相对路径。因为，永远没有绝对（）。除了根目录标记（/）和用户目录标记（～）外，有两个特定的符号用于相对路径（一个点.和两个点..）。一个点（.）代表当前工作路径，也就是用pwd执行后看到的路径；两个点（..）则代表相对当前路径的上一级目录。 比如 cd ../sub表示切换路径到上一目录的子目录 sub下；cd ../../表示切换到上一目录的上一目录；cd ./sub切换到当前目录下的子目录 sub，当然一般直接敲cd sub即可。 运行可执行文件（程序）prog：prog or ./path/prog。假如prog可被终端识别，即prog已被加入到了系统搜索路径（一般是PATH环境变量），只需要敲可执行文件名即可（比如前述cd就是系统命令，已在系统默认搜索路径里）。否则，需要加prog的路径（这里./表示执行该文件）。

强制结束终端程序运行：CTRL+C

$> sudo apt−get install python−numpy
$> pwd
$> cd /home/
$> sudodem3d

[Python 2.7]

Python与C++混合编程是目前计算软件的主流趋势。其中，Python主要扮演接口调用作用，实现快速建模；而繁重的核心计算部分则由底层C++编写的计算模块实现。SudoDEM目前仍然使用Python 2.7版本，后续项目开发升级会使用Python 3.0+版本。需要指出的是，Python 3是Python 2的全面升级，并不向下兼容，但在SudoDEM中，我们只需要最基本的功能（基本语法依然一致）。以下Python脚本展示了最为基本的入门知识。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import math # import the module math which includes some basic math functions

# we have numbers
a = 2 # integer number
b = 5
c = 2.0 # float number
d = a/b # the value is an integer
e = c/b # float number
print a, d, e # print the values of a, d and e to the terminal

#we have strings
a = 'abc' # we can assign any type of values to any type of variables
b = "this is a 'Python' script."
print a, b

#we have list, tuple and dict to store and manage the above basic primitives
c = [1,2,3,'ss'] # a list
d = (1,2,3,'ss') # a tuple
c[0] = 5 # change the first item in the list c, and the index starts from 0.
d[0] = 5 # failed! a tuple is constant, which is the distinguishing property from a list
d={1:22,2:33,4:'ssss','w':123} #a dict: 1, 2, 4 and 'w' are keys of the dict, and 22, 33, 'ssss' and 123 are the corresponding values.
print d[1], d[4], d['w']

# we initialize two objects c and d for storing data at a for loop later
c = list() # an empty list
d = dict() # an empty dict
#test expression
if 1>2:
    print "1 is greater than 2, really?"
else:
    if 2 in [1, 2, 3]: # if 2 is in the list [1, 2, 3]
        #ok, let's start a for loop
        for i in range(5): # equal to for i in [0, 1, 2, 3, 4]
            print i # you can see what's going on in the for loop
            # we append some data to a list c and a dict d
            c.append(math.cos(i)*10 + 2.0**5) # we append the value of 10*cos(i) + 2.0^5 to the list c. We use the math function cos from the module math.
            if i not in d.keys(): # if i is not a key or index of the dict d. 
                d[i] = str(i)

#we can capsule our commands by a function
def GoodJob(input_number, keyword1 = 1, keyword2 = True):
    if keyword2: # that is if keyword2 == True
        print "the input number is", input
    input_number += keyword1 # i.e., the value of input_number is updated to input_number + keyword1
    return input_number # we can return the value
    
# call the function
a = GoodJob(2) # a = 3, with print info 'the input number is 3'
b = GoodJob(3, keyword1 = 5, keyword2 = False) # b = 8, without print info.

安装SudoDEM

[二进制安装包]

在SudoDEM项目网站(https://sudodem.github.io)的Download页面下载最新的二进制安装包。或者直接到SudoDEM的Github仓库（https://github.com/SudoDEM/SudoDEM/releases）下载编译的二进制包。

将二进制包解压到本地目录（如用户目录/home/xxx/），得到以下文件目录（SudoDEM3D为例）：

通过用户配置文件.bashrc（/home/xxx/.bashrc），将可执行文件sudodem3d的路径（/home/xxx/SudoDEM/bin）加入到环境变量PATH:

PATH=${PATH}:/home/xxx/SudoDEM/bin

[脚本自动编译安装]

Here is a video that may guide you to compile SudoDEM. You can copy the command lines directly from this video.

非root权限编译安装是SudoDEM推荐安装方式。与Windows显著不同的是，Linux下向系统安装软件，会把软件的各个部分安装到系统的对应目录。比如编译的二进制会安装到bin目录（如/bin），动态链接库（为了开发等原因会把独立功能的模块编译成单独的二进制文件库，主程序根据需要动态调用该库文件）安装到lib（如 /lib）等。安装简便，但卸载可能出现问题。用sudo权限向系统安装软件（一般通过sudo apt install package），往往会附加安装一些依赖库（已存在略过）；在卸载软件时，如果选择同时卸载依赖库，那么可能导致其他软件不能正常工作（需要的依赖库被卸载）。选择非root权限安装，可以有许多好处：避免上述卸载问题、可以编译和使用特定版本的依赖库（不使用系统的）、系统更加轻量。 最新SudoDEM代码中提供了一键编译安装脚本（shell脚本），直接在终端执行即可自动下载第三方依赖库的源代码进行编译安装。首先，从SudoDEM代码仓库下载源代码，文件结构如下：

其中，文件INSTALL提供了详细安装说明，包括一键编译安装和分布手动编译安装。然后，新建第三方依赖库编译目录 3rdlib、SudoDEM2D和SudoDEM3D的编译目录build2d和build3d（按需）、SudoDEM的本地安装目录sudodeminstall，最终文件目录结构如下：

编译安装操作在scripts目录下进行，该目录包含了安装脚本如下：

第一步，安装第三方库：

./install3rdlibs.sh

如果当前系统没有安装必要的编译工具，该脚本会自动安装，并需要输入管理员密码。第三方库包括Boost、Eigen、libQGLViewer和Minieigen，SudoDEM2D和3D共用，只需编译安装一次。 第二步，安装SudoDEM程序（以SudoDEM3D为例）：

./firstCompile_SudoDEM3D.sh

注意：该脚本只用于首次编译安装；如果更改了SudoDEM的代码，只需到编译目录(build2d或build3d)执行make install即可重新编译安装。如果在SudoDEM源代码文件目录中增加了cpp文件，则需要重新cmake （即cmake ../SudoDEM3D）。

3.4 - 开源非球离散元程序SudoDEM：一图预览

颗粒形状对颗粒介质物理力学特性有着重要影响，因此，离散元非球颗粒建模有着重要意义。开源非球离散元程序SudoDEM，旨在为研究者进行非球颗粒建模提供一个快速工具。SudoDEM提供二维和三维版本，三种通用非球颗粒接触检测算法，颗粒形状涵盖球及其clump、超椭球、扩展超椭球、多面体等，Python与C++混合编程，支持OpenMP多核并行。更多功能将在后续开发中完善，欢迎大家下载、使用和测试。开发小组期待志同道合的你加入！

• 开源离散元程序SudoDEM，在第八届国际离散元大会（DEM8，荷兰，2019.7）正式对外开源发布；
• 项目主页：https://sudodem.github.io，源代码托管在GitHub仓库（https://github.com/SudoDEM）；
• 任何人在GPL v3开源协议下均可修改编辑源代码以及项目网站。

3.5 - 非球颗粒离散单元模型构建

非球颗粒离散单元模型

与球形颗粒相比，非球颗粒运动的牛顿 — 欧拉方程则采用更为一般的形式，如下： \begin{equation}\label{eqnewton_nonsph} m\frac{dv_i}{dt} = F_i \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模型则涉及组成球间的接触力计算，即颗粒内部力可求解，该模型往往被用于考虑颗粒破碎的模拟。

由球簇颗粒的定义可知，目标非球颗粒形状的逼近程度直接与组成球的个数相关。一般地，组成球的个数越多，球簇颗粒越能够逼近目标颗粒。此外，球簇颗粒中组成球的排布可以是任意的，因此，球簇模型具有很大的灵活性。理论上，球簇颗粒模型可以用于构建任意复杂形状的颗粒。但是，在实践过程中，组成球数量较大会显著地增加模拟计算耗时。鉴于此，有学者[@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}

超椭球主曲率 {#nonsph_supcurvature}

给定超椭球的参数方程为： \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。