混合有限元离散元方法

可变形离散元的坐标变换

0 接触检测算法

对非连续问题来说，判断是否接触是很重要的。减少相隔较远的颗粒间的接触判断可以提高算法效率()()o N o N →2。直接检测的算法分为两类：一类是圆形边界，另一类是矩形边界。对矩形边界，有离散单元的直接映射检测算法，单元每改变一次位置该方法需调用一次：

Step 1：将离散元映射到胞元。根据每个单元中心位置将单元映射到胞元上。因此每个单元都唯一的与一个胞元对应。即对每个单元中心坐标取整，其中d 为胞元边长，max()elements d R >

min min

int 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 设置空为-1

Loop 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 set

min

int int(

), int int()

i i i i x x i x d

y y j y d

-==+-==+11

Place the discrete element onto the corresponding singly connected list by

setting 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 set

min

int int(

), int int()

i i i i x x i x d y y j y d -==+-==+11

if C[i][j]≤N

{ C[i][j]=C[i][j]+N

loop 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 上，其中

min

int(

), int()

x y x x i d

y 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 coordinate

min

int()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 ” lists

2. 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)-list

check 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

是节点坐标，表征节点的初始位置，空间位置固定且不随节点移动，与离散单元上有限单元的边一致。基向量一般不正交，长度与有限单元相应的边的长度相同；初始形变坐标(,,)i j k

定义在节点上并随节点运动，用于表征节点初始位置向量的变形。全局、局部和初始坐标是惯性坐标，局部和初始变形坐标是非正交非单位的非惯性坐标。

局部向量和局部变形向量可以用变形梯度表示，在全局坐标下，变形张量可以写为：

(,,)(,,)(,,)u u u x

y z x u x y z v v v F f y v x y z x y z z w x y z w w w x

y z ?????+??????

?+?????????=?=?+=+???????????+????

???+???????

111 （5）在局部坐标下有

, , i i j k j i j k k i j k ??????

??????=++==++==++=??????????????????

100100001010010001 （6）

u u u u z x y x v v v v i Fi x y z x w w w w x x y z u v w i j k

x x

x ????????++??????????????????????????==+=????????????????????????????+

????????????

?????=+++ ??????11110011 … （7）通过p 点周围单元的变形f ，将局部基向量(,,)i j k 映射到局部变形向量 (,,)i j k 。变形分为()u p const =和()F p const =两种情况。()u p const =表示节点的相对位置改变，但是单元的形状体积不发生改变，这种变形不产生新的应变。

()F p const =的情况称为均匀形变，变形映射可表示为：

()()()f p f q F p q =+- （8）

例如在q 点拉伸： ()()f p q U p q =+- （9）其中U 是对称正定方阵，设U 的张量基为(,,)e e e 123，则U 可以表示成对角阵

s U s s ??

??=?

????

3 （11）其中(,,)s s s 123是特征值，一般在不同的方向拉伸或缩短的长度是不同的。

在q 点旋转： ()()f p q R p q =+- （12）

其中R 是单位正交阵，即 det()T R R R -==11且。对任意向量a

()R Ra R Ra Ia --==11

记转轴为ψ，旋转角为ψ，则任取一组含ψ的正交基(,,

)e e ψ

12，R 可以表示为 cos sin sin cos R ψ

ψψ

-??

??=?????

1 (13)

2 均匀形变

综上（8）—（13）我们可将均匀形变（F=const ）表示为旋转与拉伸的组合：

()f p g s s g F RU VR

==== （14）

其中梯度R ，U ，V 都是常数张量，不随节点的不同而改变。

det det det U V F ==>0 （15）

U 是右乘拉伸张量，表示先做拉伸后做旋转，可表示为三个坐标方向的拉伸叠加：

i i , U (), V ()i i i i i i i i i i i i i

i U U U U e e I s e e V V V V e e I s e e λλ====?=+-?==?=+-?∑∑3

12313

123111

其中i e 是局部基向量， i e 是局部变形基向量， i i e Re =。应变：可定义不同的应变张量，如Cauchy-Green 应变张量

T T T T C F F U R RU U U ===

设伸长量为s ，则变形量为

n ln() n n s e n

s ?-≠?

=??=?

00 （16）应变张量的定义格式随n 取值而不同。

应力：由柯西应力的定义对向量m ，表面张力S(m)=Tm 。T 是对称正定的，且满足运动方程，

divT b v

ρ+= （17）其中S 是变形区域产生的张力，ρ单位体积变形区域的密度，向量m 的方向垂直于变形表面，长度等于变形域的面积。对于全局坐标，

xy xz yx

yy yz zx zy

zz t t t T t t t t t t ????=?????

()x x xx

xy xz x y z y yx

yy yz y zx zy

zz z z S m t t t S m S i S j S k S t t t m t t t S m ????

??????

?=++==?????????????

?????

（18）

3 单元变形（应变为常数）以下考虑二维三角形单元，

初始变形坐标可以用初始坐标表示：

, x y x y i i i i j j j i j j =+=+

单元上任意向量可以通过坐标变换转到不同的坐标系下。三角单元的变形由节点处的变形表示，因为F=const ，所以单元变形可用线性函数表示，令(,)i i x y 是全局坐标下变形前的分量，(,)f f x y 是变形后的，则

, f x i x i f y i y i x x y y x y αβαβ=+=+ （19）

4 中心差分时间积分算法

如下是显式时间积分算法及迭代格式： Verlet 算法（x 是位移，v 是速度）： h t =?

, t h t h t t h t t h h

x x hv v v f m

+-++=+=+

22 T-1/6算法：

t t h t t t t t t t t t h t h f x x v h h b v h m f v h b v h m v b h +++=++

+++=

-232

261

其中() and ()t t t h

t h df df b b m dx m dx

++==11。此外还有T-1/12，D-1/12算法。上述方法都有相同的迭代格式，由三步构成： (1) 预测：由t 时刻的位移及其导数计算t+h 时刻的位移

,,,!!!

!!...!

iv t h p t t t t t iv t h p t t t t iv t h p

t t t h h h x x x h x x x h h x x x h x x h x x x h x +++=++++=+++=++234

2234232 (2) 估计：在这一阶段，在t+h 时刻的力t h f +由预测步中的,t h p x +值计算得到。然

后由t h f +求得t+h 时刻的加速度t h x +

，计算差值 ,t h t h p x x x

++?=- (3) 校正：校正位移

,,,,!!

!!!...

!!!

t h t h p t h t h p t h t h p t h t h p xh x x xh x

h x h h h xh x x h h xh x x αααα++++++++?=+?=+?=+?=+2

1222

332

322222332 对三阶格式，/, /, , /αααα====01231656113。