基于DIC的颗粒间接触力计算及力链分析_陈凡秀
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
4.2 颗粒平均应力的计算
由于单层钢珠的厚度相对于整个平面面积可 以忽略，加载外力及容器壁约束力均在钢珠平面内， 因此本次试验的应力分析可简化为平面应力问题。 根据平面应力问题的物理方程(8)，该方程可推导为 式(9)。
均应力方程。以 Kα 矩阵对应的计算为例，前两行 为力的平衡，第三行为力矩平衡，最后三行为平均 应力方程。对矩阵方程组求解即可求得颗粒间每个 接触点对应的接触力的大小和方向。
246
应用力学学报
第 32 卷
展开如下
α
β
⎡% 0 " 0 "⎤ ⎧ # ⎫ ⎧ #
⎫
p
⎢ ⎢
0
⎢#
Kα 0
0 Kβ %0
0 #
⎥ ⎥ ⎥
⎪ ⎪⎪ ⎨
fα #
⎪ ⎪⎪ ⎬
=
⎪⎪⎪b p ⎨
−K #
β
f
β
⎪ ⎪⎪ ⎬
q⎢⎢" −K α "
0
"⎥⎥
⎪ ⎪
f
β
⎪ ⎪
⎪ ⎪
bq
⎪ ⎪
⎢⎣ # 0 # 0 %⎥⎦ ⎪⎩ # ⎪⎭ ⎪⎩ #
移，并获得εx、εy、γxy。对图像进行边缘提取得到圆 心坐标和半径，对每个颗粒面积区域内的应力求和
4 二维钢珠颗粒受集中力试验
取平均，得到各颗粒的平均应力。
4.1 试验系统
二维钢珠颗粒受集中力的试验系统见图2。 CCD相机采用的是德国Basler产品(scA1600-14fm)。 施加压力的装置为岛津牌万能试验机，型号为
图 4 钢珠排列结构 Fig.4 Information on the fabric: particle number, position of
到真实颗粒中接触力的大小及形成力链时接触力的 变化规律，直观地观察颗粒间力的传递及力链的演 变过程。
2 数字图像相关方法
DIC[12]是一种基于计算机视觉原理、图像处理、 数值运算的一种非接触、非干涉的光测力学方法， 该方法中的设备简单且抗干扰能力强，在当前实验 力学领域具有广泛的应用。它的基本原理是对变形 前后的图像子区进行相关性计算，图像子区一般选 取以待求点为中心点的正方形区域。对变形后的子 区按恰当的相关函数进行计算，选取相关性最大的 变形后子区为目标图像子区，来计算图像的全场应 变。选用标准化协方差相关函数[13]，该式为
1引言
在工程实际中，研究颗粒物质的力学性质具 有非常重要的意义。颗粒间力的传播路径（即力 链），决定着颗粒物质宏观及微观的力学性质。近 年来，研究人员通过各种方法对颗粒受力进行研 究。文献[1]研究了二维光弹颗粒体系力链长度的 分布，并得到力链长度随施加力的变化呈指数衰 减的规律。文献[2]使用光弹性测量技术将二维颗粒 中的受力及力的传递过程可视化并对力的大小进行 定量测量。文献[3]用离散元模拟对二维颗粒材料的 受力变形进行了研究，得到内部的力链屈曲引起宏 观上的剪胀现象的结论。文献[4]通过模拟实验研究
为接触点到原点的坐标矢量。如图 1 所示，与 P 颗
粒接触的 Ncp 个接触力在颗粒 P 的面积 Ωp 范围内满
足力的平衡，这些接触力包括颗粒与颗粒之间的接
触力 f α 、已知的边界力 f β 和未知的边界力 f γ 。
xα 表示第 α 个接触点到图像坐标原点的坐标矢量， 力矩平衡即 P 颗粒上的所有接触力到图像原点的力 矩满足平衡。
245
料的不透明材料颗粒间的接触力。 通过牛顿力学理论和颗粒线性动量平衡求颗
粒间接触力的方法，在真实颗粒试验中的应用研究 较少，本文以二维钢珠受集中压力的试验为例，对 真实颗粒的受力进行研究。在二维钢珠颗粒体系受 钢条施加的集中压力时，采用 CCD 相机连续拍摄 二 维 图 像 ， 通 过 数 字 图 像 相 关 方 法 [8-10] (Digital Image Correlation， DIC)得到钢条在向下移动过程 中钢珠颗粒体系内颗粒受力时产生的位移和应变， 从而获得全场应力值，并求得每个颗粒的平均应力。 对采集的初始图像进行边缘提取及霍夫(Hough)变 换[11]，得到每个颗粒的形心位置坐标和半径以及颗 粒间接触点的位置坐标；对所有颗粒和接触点进行 编号，获得颗粒体系的排列结构。将平均应力和颗 粒体系排列结构作为初始数据，通过牛顿力学理论 和颗粒线性动量平衡求得颗粒间接触力，定量地得
基金项目：国家自然科学基金(11472145；51008166)
收稿日期：2014-10-30
修回日期：2015-03-11
第一作者简介：陈凡秀，女，1979 年生，博士，青岛理工大学，副教授；研究方向——光测力学的理论与应用。 E-mail：mecfx@163.com
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
AG-IC50kN。该试验机可实时显示钢条底部在钢条 下压过程中的受力并通过相连的计算机绘制成钢条 底部接触力随压入深度的变化图。试验采用直径为 5mm的钢珠，为了防止钢珠表面反光并区分灰度， 在钢珠表面人工形成散斑。在凹形槽中铺入单层的 钢珠颗粒，用透明有机玻璃覆盖并将有机玻璃固定 于凹形槽上作为观测面。将与万能试验机压力探头 连接的钢条(横截面积：1mm×9mm)以恒定速率压 入钢珠颗粒体系中，CCD相机以2幅/秒的频率采集 体系变形图。
4.3 颗粒体系的排列结构
由于钢珠上喷涂了散斑，直接对原始图片提取 钢珠边缘较难实现，需要对图像进行初步的处理。 首先对初始图像进行二值化处理，如图 3(b)所示， 目的是将钢珠颗粒与背景更好地区分开。采用自适
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
247
应 Canny 算子[14]对二值化后的图像进行初步边缘提 取，Canny 算子具有较好的抗噪性能和低误码率， 提取结果见图 3(c)。对初步提取出钢珠边缘的图像 采用 Hough 圆变换进行检测。Hough 圆变换是将平
文章编号：1000- 4939(2015) 02-0244-07
Vol.32 No.2 Apr. 2015
基于 DIC 的颗粒间接触力计算及力链分析
陈凡秀 张慧新 庄琦
(青岛理工大学 理学院 266033 青岛)
摘要：首次利用数字图像相关方法(Digital Image Correlation method，DIC)，依据牛顿力学理论和 颗粒线性动量平衡，分析了集中力作用下的二维钢珠颗粒体系，定量地获得了颗粒间接触力。采 用 CCD 相机进行观测并采集了集中力作用下钢珠颗粒体系的变形图像；利用数字图像相关方法 对采集的序列图像进行处理，获得了全场位移及应变，并根据胡克定律求得全场应力。对每个颗 粒面积区域内各点应力值求和并取平均，获得各个颗粒的平均应力。对采集的数字图像进行图像 分析，得到各个颗粒及接触点的排列结构。以颗粒平均应力和颗粒的排列结构作为初始条件，根 据牛顿力学理论和颗粒线性动量平衡，最终获取了颗粒间接触力的大小和方向，实现了真实颗粒 体系接触力的定量计算。另外从受力方面对集中力作用下颗粒体系内部力链的发展与演变进行了 分析。研究结果表明：力链起源于加载位置，随着外荷载的部分耗散，颗粒间接触力逐渐发生变 化，力链向四周扩展而形成力链网络。通过对不同压入深度颗粒间接触力变化和分布的分析，定 量地描述了加载过程中颗粒间的滚动或滑动引起的力链断裂和重构的现象。 关键词：颗粒体系；数字图像相关方法；受力分析；接触力；力链 中图分类号：TU432；TU13 文献标识码：A DOI: 10.11776/cjam.32.02.D055
P 颗粒上的平均应力为对 Ωp 面积范围内的应
力值求平均，即
图 1 颗粒边界力和颗粒间接触力 Fig.1 Illustration of particle-particle and particle-boundary contacts
考虑线性动量平衡、发散定理和柯西应力张量
的对称性，公式(4)可写为
的圆心和接触点分别编号，得到颗粒整体的排列结
构，如图 4 所示。其中粗体为圆心编号，细体为接
触点编号，图中共有 186 个颗粒和 435 个接触点。 通过所有的圆心和接触点形成颗粒排列结构，以获 取坐标矢量 xα 和各颗粒的平均应力。将得到的平均应 力值和颗粒排列结构作为原始数据，根据式(7)即可求 得颗粒间接触力。图 5 给出了接触力计算的流程图。
面圆边缘上的点对应到参数空间中相应的圆上，在
参数空间内对点进行描述，对可能会落在边缘上的
点进行统计计算，根据统计数据的结果确定属于边
缘的程度，精确检测出每个钢珠颗粒的圆心坐标和
半径，见图 3(d)。还需要根据检测出的圆心坐标和 半径获得颗粒整体的排列结构。以两圆心距离在颗
粒直径范围内为限定条件求出接触点坐标。对所有
网络出版时间：2015-04-20 11:43 网络出版地址：http://www.cnki.net/kcms/detail/61.1112.O3.20150420.1143.026.html
第 32 卷 第 2 期 2015 年 4 月
应用力学学报 CHINESE JOURNAL OF APPLIED MECHANICS
σx、σy、τxy 分别为 x 方向、y 方向、切向应力；E 为
弹性模量；μ 为泊松比。
根据颗粒大小和喷制散斑的情况将图像子区
大小设为47pixel，为求得全场每个像素点的位移将 步长设为1pixel。用数字图像处理软件对钢条下压过 程中采集的一系列图像进行相关处理，得到全场位
图 2 试验系统 Fig.2 Experimental system
MM
∑ ∑ [ f (x, y) − fm ][g(x + u, y + v) − gm ]
C(u,v) =
x=−M y=−M
(1)
MM
MM
∑ ∑ [ f (x, y) − ]fm 2 ∑ ∑ [g(x + u, y + v) − gm ]2
x=−M y=−M
x=−M y=−M
式中：f(x,y)、g(x+u,y+v)分别为变形前后数字图像 中各像素点的灰度值；fm、gm 为其图像子区的平均 灰度值；u、v 为子区中心的位移，以像素为单位。
∑ ( ) σ p =
1
Ncp
sym
f α ⊗ xα
Ωp α =1
(5)
其中： f α 为待求的颗粒间接触力；符号 ⊗ 表示并
列运算，比如 (α ⊗ β )ijkl = αij βkl 。
在二维平面上，式(2)、式(3)、式(5)可用矩阵形 式表示，即
Ncp
∑ K α ⋅ f α =b p
(6)
α =1
了颗粒物质受力时的性质，提出颗粒物质具有多尺 度特性，并对颗粒物质宏观、细观、微观尺度上的 规律做了统计分析；提出颗粒间判断强力链的大小 和角度方法[5]，分析了压入试验过程中力链的演变 规律。
对有缺陷和不确定性的真实颗粒内部力的大 小、分布、力链提取及分析尚未有合适的研究方法。 2012 年，Andrade 和 Avila[6]提出了 GEM(Granular Element Method)方法，GEM 是基于牛顿力学理论和 颗粒线性动量平衡而建立的一种求颗粒间接触力的 方法。文中通过非接触的 FEM 方法得到初始数据， 包括颗粒排列结构、宏观集中力与颗粒本身所受的 应力，经过 GEM 计算得到细观的颗粒间接触力。 Ryan Hurley[7]等利用 GEM 方法研究了喷涂光弹涂
⎪⎭
(7) 式中
⎡1
⎢⎢0
K
α
=
⎢ ⎢ ⎢
x2α x1α
⎢
⎢0
⎢⎣ x2α
0⎤
1
⎥ ⎥
− x1α 0 x2α
⎥ ⎥ ⎥ ⎥ ⎥
，
f
α
=
⎧⎪ ⎨
⎪⎩
f
α
1
f
α
2
⎫⎪ ⎬
，
⎪⎭
x1α ⎥⎦
⎧0⎫
⎪ห้องสมุดไป่ตู้⎪
0
⎪ ⎪
⎪0⎪
bp
=
⎪ ⎨⎪Ωp
σ
p 11
⎪ ⎬ ⎪
⎪⎪Ωp
σ
p 22
⎪ ⎪
⎪⎩2Ωp
σ
p 12
⎪ ⎭
对每个颗粒都可建立力的平衡、力矩平衡、平
3 颗粒间接触力计算
∫ σ p = 1
Ωp
Ωp σ dΩp
(4)
其中：σ p 为平均应力；σ 为每个像素点上的应力。
根据牛顿力学定律确定力的平衡和力矩平衡，
见式(2)和式(3)，即
Ncp
∑ f α =0
(2)
α =1
Ncp
∑ f α × xα =0
(3)
α =1
其中：Ncp 为接触点个数； f α 为颗粒间接触力；xα
⎧⎪ε ⎪
x
=
1（σ E
x
-μσ
）
y
⎪⎨ε ⎪
y
=
1（σ E
y
-μσ
）
x
(8)
⎪ ⎪⎩
γ
xy
=
（2 1+μ）τ E
xy
⎧⎪σ ⎪
x
（ε =
x +με y）E 1− μ2
⎪⎪⎨σ ⎪
y
=（με1x−+εμy2）E
(9)
⎪ ⎪ ⎪⎩
τ
xy
=
（2 1E+μ）γ
xy
其中：εx、εy、γxy 分别为 x 方向、y 方向、切向应变；
由于单层钢珠的厚度相对于整个平面面积可 以忽略，加载外力及容器壁约束力均在钢珠平面内， 因此本次试验的应力分析可简化为平面应力问题。 根据平面应力问题的物理方程(8)，该方程可推导为 式(9)。
均应力方程。以 Kα 矩阵对应的计算为例，前两行 为力的平衡，第三行为力矩平衡，最后三行为平均 应力方程。对矩阵方程组求解即可求得颗粒间每个 接触点对应的接触力的大小和方向。
246
应用力学学报
第 32 卷
展开如下
α
β
⎡% 0 " 0 "⎤ ⎧ # ⎫ ⎧ #
⎫
p
⎢ ⎢
0
⎢#
Kα 0
0 Kβ %0
0 #
⎥ ⎥ ⎥
⎪ ⎪⎪ ⎨
fα #
⎪ ⎪⎪ ⎬
=
⎪⎪⎪b p ⎨
−K #
β
f
β
⎪ ⎪⎪ ⎬
q⎢⎢" −K α "
0
"⎥⎥
⎪ ⎪
f
β
⎪ ⎪
⎪ ⎪
bq
⎪ ⎪
⎢⎣ # 0 # 0 %⎥⎦ ⎪⎩ # ⎪⎭ ⎪⎩ #
移，并获得εx、εy、γxy。对图像进行边缘提取得到圆 心坐标和半径，对每个颗粒面积区域内的应力求和
4 二维钢珠颗粒受集中力试验
取平均，得到各颗粒的平均应力。
4.1 试验系统
二维钢珠颗粒受集中力的试验系统见图2。 CCD相机采用的是德国Basler产品(scA1600-14fm)。 施加压力的装置为岛津牌万能试验机，型号为
图 4 钢珠排列结构 Fig.4 Information on the fabric: particle number, position of
到真实颗粒中接触力的大小及形成力链时接触力的 变化规律，直观地观察颗粒间力的传递及力链的演 变过程。
2 数字图像相关方法
DIC[12]是一种基于计算机视觉原理、图像处理、 数值运算的一种非接触、非干涉的光测力学方法， 该方法中的设备简单且抗干扰能力强，在当前实验 力学领域具有广泛的应用。它的基本原理是对变形 前后的图像子区进行相关性计算，图像子区一般选 取以待求点为中心点的正方形区域。对变形后的子 区按恰当的相关函数进行计算，选取相关性最大的 变形后子区为目标图像子区，来计算图像的全场应 变。选用标准化协方差相关函数[13]，该式为
1引言
在工程实际中，研究颗粒物质的力学性质具 有非常重要的意义。颗粒间力的传播路径（即力 链），决定着颗粒物质宏观及微观的力学性质。近 年来，研究人员通过各种方法对颗粒受力进行研 究。文献[1]研究了二维光弹颗粒体系力链长度的 分布，并得到力链长度随施加力的变化呈指数衰 减的规律。文献[2]使用光弹性测量技术将二维颗粒 中的受力及力的传递过程可视化并对力的大小进行 定量测量。文献[3]用离散元模拟对二维颗粒材料的 受力变形进行了研究，得到内部的力链屈曲引起宏 观上的剪胀现象的结论。文献[4]通过模拟实验研究
为接触点到原点的坐标矢量。如图 1 所示，与 P 颗
粒接触的 Ncp 个接触力在颗粒 P 的面积 Ωp 范围内满
足力的平衡，这些接触力包括颗粒与颗粒之间的接
触力 f α 、已知的边界力 f β 和未知的边界力 f γ 。
xα 表示第 α 个接触点到图像坐标原点的坐标矢量， 力矩平衡即 P 颗粒上的所有接触力到图像原点的力 矩满足平衡。
245
料的不透明材料颗粒间的接触力。 通过牛顿力学理论和颗粒线性动量平衡求颗
粒间接触力的方法，在真实颗粒试验中的应用研究 较少，本文以二维钢珠受集中压力的试验为例，对 真实颗粒的受力进行研究。在二维钢珠颗粒体系受 钢条施加的集中压力时，采用 CCD 相机连续拍摄 二 维 图 像 ， 通 过 数 字 图 像 相 关 方 法 [8-10] (Digital Image Correlation， DIC)得到钢条在向下移动过程 中钢珠颗粒体系内颗粒受力时产生的位移和应变， 从而获得全场应力值，并求得每个颗粒的平均应力。 对采集的初始图像进行边缘提取及霍夫(Hough)变 换[11]，得到每个颗粒的形心位置坐标和半径以及颗 粒间接触点的位置坐标；对所有颗粒和接触点进行 编号，获得颗粒体系的排列结构。将平均应力和颗 粒体系排列结构作为初始数据，通过牛顿力学理论 和颗粒线性动量平衡求得颗粒间接触力，定量地得
基金项目：国家自然科学基金(11472145；51008166)
收稿日期：2014-10-30
修回日期：2015-03-11
第一作者简介：陈凡秀，女，1979 年生，博士，青岛理工大学，副教授；研究方向——光测力学的理论与应用。 E-mail：mecfx@163.com
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
AG-IC50kN。该试验机可实时显示钢条底部在钢条 下压过程中的受力并通过相连的计算机绘制成钢条 底部接触力随压入深度的变化图。试验采用直径为 5mm的钢珠，为了防止钢珠表面反光并区分灰度， 在钢珠表面人工形成散斑。在凹形槽中铺入单层的 钢珠颗粒，用透明有机玻璃覆盖并将有机玻璃固定 于凹形槽上作为观测面。将与万能试验机压力探头 连接的钢条(横截面积：1mm×9mm)以恒定速率压 入钢珠颗粒体系中，CCD相机以2幅/秒的频率采集 体系变形图。
4.3 颗粒体系的排列结构
由于钢珠上喷涂了散斑，直接对原始图片提取 钢珠边缘较难实现，需要对图像进行初步的处理。 首先对初始图像进行二值化处理，如图 3(b)所示， 目的是将钢珠颗粒与背景更好地区分开。采用自适
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
247
应 Canny 算子[14]对二值化后的图像进行初步边缘提 取，Canny 算子具有较好的抗噪性能和低误码率， 提取结果见图 3(c)。对初步提取出钢珠边缘的图像 采用 Hough 圆变换进行检测。Hough 圆变换是将平
文章编号：1000- 4939(2015) 02-0244-07
Vol.32 No.2 Apr. 2015
基于 DIC 的颗粒间接触力计算及力链分析
陈凡秀 张慧新 庄琦
(青岛理工大学 理学院 266033 青岛)
摘要：首次利用数字图像相关方法(Digital Image Correlation method，DIC)，依据牛顿力学理论和 颗粒线性动量平衡，分析了集中力作用下的二维钢珠颗粒体系，定量地获得了颗粒间接触力。采 用 CCD 相机进行观测并采集了集中力作用下钢珠颗粒体系的变形图像；利用数字图像相关方法 对采集的序列图像进行处理，获得了全场位移及应变，并根据胡克定律求得全场应力。对每个颗 粒面积区域内各点应力值求和并取平均，获得各个颗粒的平均应力。对采集的数字图像进行图像 分析，得到各个颗粒及接触点的排列结构。以颗粒平均应力和颗粒的排列结构作为初始条件，根 据牛顿力学理论和颗粒线性动量平衡，最终获取了颗粒间接触力的大小和方向，实现了真实颗粒 体系接触力的定量计算。另外从受力方面对集中力作用下颗粒体系内部力链的发展与演变进行了 分析。研究结果表明：力链起源于加载位置，随着外荷载的部分耗散，颗粒间接触力逐渐发生变 化，力链向四周扩展而形成力链网络。通过对不同压入深度颗粒间接触力变化和分布的分析，定 量地描述了加载过程中颗粒间的滚动或滑动引起的力链断裂和重构的现象。 关键词：颗粒体系；数字图像相关方法；受力分析；接触力；力链 中图分类号：TU432；TU13 文献标识码：A DOI: 10.11776/cjam.32.02.D055
P 颗粒上的平均应力为对 Ωp 面积范围内的应
力值求平均，即
图 1 颗粒边界力和颗粒间接触力 Fig.1 Illustration of particle-particle and particle-boundary contacts
考虑线性动量平衡、发散定理和柯西应力张量
的对称性，公式(4)可写为
的圆心和接触点分别编号，得到颗粒整体的排列结
构，如图 4 所示。其中粗体为圆心编号，细体为接
触点编号，图中共有 186 个颗粒和 435 个接触点。 通过所有的圆心和接触点形成颗粒排列结构，以获 取坐标矢量 xα 和各颗粒的平均应力。将得到的平均应 力值和颗粒排列结构作为原始数据，根据式(7)即可求 得颗粒间接触力。图 5 给出了接触力计算的流程图。
面圆边缘上的点对应到参数空间中相应的圆上，在
参数空间内对点进行描述，对可能会落在边缘上的
点进行统计计算，根据统计数据的结果确定属于边
缘的程度，精确检测出每个钢珠颗粒的圆心坐标和
半径，见图 3(d)。还需要根据检测出的圆心坐标和 半径获得颗粒整体的排列结构。以两圆心距离在颗
粒直径范围内为限定条件求出接触点坐标。对所有
网络出版时间：2015-04-20 11:43 网络出版地址：http://www.cnki.net/kcms/detail/61.1112.O3.20150420.1143.026.html
第 32 卷 第 2 期 2015 年 4 月
应用力学学报 CHINESE JOURNAL OF APPLIED MECHANICS
σx、σy、τxy 分别为 x 方向、y 方向、切向应力；E 为
弹性模量；μ 为泊松比。
根据颗粒大小和喷制散斑的情况将图像子区
大小设为47pixel，为求得全场每个像素点的位移将 步长设为1pixel。用数字图像处理软件对钢条下压过 程中采集的一系列图像进行相关处理，得到全场位
图 2 试验系统 Fig.2 Experimental system
MM
∑ ∑ [ f (x, y) − fm ][g(x + u, y + v) − gm ]
C(u,v) =
x=−M y=−M
(1)
MM
MM
∑ ∑ [ f (x, y) − ]fm 2 ∑ ∑ [g(x + u, y + v) − gm ]2
x=−M y=−M
x=−M y=−M
式中：f(x,y)、g(x+u,y+v)分别为变形前后数字图像 中各像素点的灰度值；fm、gm 为其图像子区的平均 灰度值；u、v 为子区中心的位移，以像素为单位。
∑ ( ) σ p =
1
Ncp
sym
f α ⊗ xα
Ωp α =1
(5)
其中： f α 为待求的颗粒间接触力；符号 ⊗ 表示并
列运算，比如 (α ⊗ β )ijkl = αij βkl 。
在二维平面上，式(2)、式(3)、式(5)可用矩阵形 式表示，即
Ncp
∑ K α ⋅ f α =b p
(6)
α =1
了颗粒物质受力时的性质，提出颗粒物质具有多尺 度特性，并对颗粒物质宏观、细观、微观尺度上的 规律做了统计分析；提出颗粒间判断强力链的大小 和角度方法[5]，分析了压入试验过程中力链的演变 规律。
对有缺陷和不确定性的真实颗粒内部力的大 小、分布、力链提取及分析尚未有合适的研究方法。 2012 年，Andrade 和 Avila[6]提出了 GEM(Granular Element Method)方法，GEM 是基于牛顿力学理论和 颗粒线性动量平衡而建立的一种求颗粒间接触力的 方法。文中通过非接触的 FEM 方法得到初始数据， 包括颗粒排列结构、宏观集中力与颗粒本身所受的 应力，经过 GEM 计算得到细观的颗粒间接触力。 Ryan Hurley[7]等利用 GEM 方法研究了喷涂光弹涂
⎪⎭
(7) 式中
⎡1
⎢⎢0
K
α
=
⎢ ⎢ ⎢
x2α x1α
⎢
⎢0
⎢⎣ x2α
0⎤
1
⎥ ⎥
− x1α 0 x2α
⎥ ⎥ ⎥ ⎥ ⎥
，
f
α
=
⎧⎪ ⎨
⎪⎩
f
α
1
f
α
2
⎫⎪ ⎬
，
⎪⎭
x1α ⎥⎦
⎧0⎫
⎪ห้องสมุดไป่ตู้⎪
0
⎪ ⎪
⎪0⎪
bp
=
⎪ ⎨⎪Ωp
σ
p 11
⎪ ⎬ ⎪
⎪⎪Ωp
σ
p 22
⎪ ⎪
⎪⎩2Ωp
σ
p 12
⎪ ⎭
对每个颗粒都可建立力的平衡、力矩平衡、平
3 颗粒间接触力计算
∫ σ p = 1
Ωp
Ωp σ dΩp
(4)
其中：σ p 为平均应力；σ 为每个像素点上的应力。
根据牛顿力学定律确定力的平衡和力矩平衡，
见式(2)和式(3)，即
Ncp
∑ f α =0
(2)
α =1
Ncp
∑ f α × xα =0
(3)
α =1
其中：Ncp 为接触点个数； f α 为颗粒间接触力；xα
⎧⎪ε ⎪
x
=
1（σ E
x
-μσ
）
y
⎪⎨ε ⎪
y
=
1（σ E
y
-μσ
）
x
(8)
⎪ ⎪⎩
γ
xy
=
（2 1+μ）τ E
xy
⎧⎪σ ⎪
x
（ε =
x +με y）E 1− μ2
⎪⎪⎨σ ⎪
y
=（με1x−+εμy2）E
(9)
⎪ ⎪ ⎪⎩
τ
xy
=
（2 1E+μ）γ
xy
其中：εx、εy、γxy 分别为 x 方向、y 方向、切向应变；