混合有限元离散元方法

可变形离散元的坐标变换0 接触检测算法对非连续问题来说，判断是否接触是很重要的。

减少相隔较远的颗粒间的接触判断可以提高算法效率()()o N o N →2。

直接检测的算法分为两类：一类是圆形边界，另一类是矩形边界。

对矩形边界，有离散单元的直接映射检测算法，单元每改变一次位置该方法需调用一次：Step 1：将离散元映射到胞元。

根据每个单元中心位置将单元映射到胞元上。

因此每个单元都唯一的与一个胞元对应。

即对每个单元中心坐标取整，其中d 为胞元边长，max()elements d R >min minint int(), int int()i i i i x x y y x y d d--=+=+11 （1） Step 2：找到可能接触的离散元。

每个离散元都与其映射胞元中心重合且被完全覆盖，两个胞元相邻的单元有可能接触。

事实上对下面的胞元，因为每个胞元都可作为中心元被检查两次，所以只需检查左下角的四个。

如何将单元映射到胞元上呢？不同的映射方法定义了不同的算法，如屏风法、分类法、Munjiza-NBS 法和Williams-C 网格法。

1 屏风法：[,,,...,]N E e e e e =123，N 为单元个数。

离散元素映射到细胞的获得是通过设置阵列C 指单独列出的头连接的胞元。

映射进行如下： STEP 1 设置空为-1Loop over all rows of cells(i=1; i ≤Ncel ; i++){ Loop over all cells in a row(j=1;j ≤Ncel ; j++) { set C[i][j ] = −1 } }Loop over all discrete elements(k=1; k ≤N; k++) { set E[k] = −1 }STEP 2 离散元映射到胞元。

单元每运动一次执行一次，例如每个时间步执行一次。

对每个胞元设置一个独立的链表：Loop over all discrete elements (k=1; k ≤N;k++) { Integerise current coordinates and setminminint int(), int int()i i i i x x i x dy y j y d-==+-==+11Place the discrete element onto the corresponding singly connected list bysetting E[k] = C[i][j ] and afterwards C[i][j ] = k}当单元移动后再进行映射时不需要初始化。

当单元分离时步骤1再次进行。

一旦 离散元都映射到胞元，则检测开始。

对所有离散元循环，找到离散元映射的胞元，然后检查该胞元相邻和周围胞元中的单元。

Loop over discrete elements (k=1; k ≤N; k++) { Integerise current coordinates and setminminint int(), int int()i i i i x x i x d y y j y d -==+-==+11if C[i][j]≤N{ C[i][j]=C[i][j]+Nloop over all discrete elements from C[i][j] list{ loop over all discrete elements from neighboring cells, i.e. lists C[i-1][j-1], C[i][j-1], C[i+1][j-1], C[i-1][j], C[i][j]{ direct check for contact between discrete elements } } }}搜索只遍历所有离散元，并不遍历所有胞元。

2 分类法如果离散元是分散分布的，则胞元矩阵会很大，所以为节省RAM 和CPU 产生了新的方法。

设置数组X 记录胞元所在的列，数组Y 记录胞元所在的列，数组D 记录离散元，长度皆为N 。

3 Munjiza-NBS 法该算法是找到相隔最近的点对，即它们最近的点之间的距离小于或等于零，则单元重叠或接触。

空间重构：空间被分为理想的长度为d 的正方形胞元。

从集合{,,...,}E N =12映射到(,)(,)(,)...(,)(,)(,)(,)...(,)...............(,)(,)(,)...(,)n n C m m m m n ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦11121312122232123 将离散元(x,y)映射到胞元(,)x y i j 上，其中minminint(), int()x y x x i dy y i d -=-=STAGE1 链接列：对每一列y i (y y i N ≤)，构造一个独立的接触链，包含在这一列上的所有离散元。

单元序号按降序排列相互关联，例→→→10741，它们对应的胞元y i 都等于7。

Loop over all discrete elements ( i=1; i ≤N; i++) { calculate integerised coordinateminint()y y y i d -=and place the discrete element onto the y yi list }数组B 存放的元素是最后一个映射进当前列链的单元，有N y 个分量，指向数组Y 的相应位置，引出一个新的链表。

数组Y 中存放着当前单元指向的单元，当指向结束时当前值设为-1，有N 个分量。

STAGE 2：对所有y y i =的单元循环，将x x i =的单元放入链表X ，且令该链表为NEW 。

有N 个分量。

检查一个链表后令其为OLD 。

则每建立一个新链表都应用旧链表的空间，只是需先进行初始化。

算法总结：1. Loop over all discrete elements( k=0; k ≤N; k++){ calculate integerised y coordinates of its centre place the discrete element onto a list for the corresponding row of cells (y-list) }consider all y-lists to be “new ” lists2. Loop over discrete elements( k=0; k ≤N; k++){ calculate integerised y coordinates of its centre if discrete element belongs to a new y-list{ mark the y-list as an old list and call it central y-list { 3. loop over all discrete elements from central y-list{ integerise x coordinate of the discrete element and place it onto the corresponding x-list }4. loop over all discrete elements from neighbouring (y-)-list{ integerise x coordinate of the discrete element and place discrete element onto the corresponding x-list }5. loop over all discrete elements from the central y-list { if the discrete element belongs to a new list (x,y) list { mark the list (x,y) as old and call it central (x,y)-listcheck for contact between discrete elements from the central list (x,y) and lists (x,y),(x-1,y),(x-1,y-1),(x,y-1), (x+1,y-1) } } }6. loop over all discrete element from the central y-list{ remove the corresponding x-list, i.e. set the head of the list to zero } 7. loop over all discrete element from list (y-1){ remove the corresponding x-list, i.e. set head of the list to zero } } }（2）8. loop over all discrete elements{ remove the corresponding y-lists i.e. set the head of the list to zero}这种方法的计算时间 O(N) ，需要的RAM 空间为A[2][x n ]，B[y n ]，X[N]和Y[N]，共x y M n n N =++22 个。

1 基本概念每个离散的单元代表一个独立的变形体，在任一时间点占据一定的空间区域。

对这些区域赋予特殊的意义，例如单元上所有点在T=0时刻的初始位置构成一个集合，定义为区域B 。

区域B 的任意子域称作该单元的一个块，单元体的变形通过如下映射定义：()x f p = (1)f 是将B 映射到闭区域E 的光滑双射，且满足()f p ∇>0，Ker f =Null (2)(1)式还可表示为 ()()x f p p u p ==+， (3)其中()u p 为位移，映射 ()()F p f p I u =∇=+∇ (4) 称为变形梯度，描述一点附近的形变量。

为了描述在离散单元上一点P 附近物质的变形，定义四组坐标： 全局坐标 (i,j,k )；局部坐标(,,)i j k 用于描述离散单元的初始位置，因此空间位置固定且不随单元移动，一般同单元的重力方向一致；局部变形坐标 (,,)i j k定义在节点上并随节点运动，表征节点附近的变形。

方向与大小都随单元的变形而改变，遵循(1)式，在初始点是两两正交的非单位向量，变形后不再保持正交性；初始坐标(,,)i j k是节点坐标，表征节点的初始位置，空间位置固定且不随节点移动，与离散单元上有限单元的边一致。