基于DIC的颗粒间接触力计算及力链分析_陈凡秀

到真实颗粒中接触力的大小及形成力链时接触力的 变化规律，直观地观察颗粒间力的传递及力链的演 变过程。
2 数字图像相关方法
DIC[12]是一种基于计算机视觉原理、图像处理、 数值运算的一种非接触、非干涉的光测力学方法， 该方法中的设备简单且抗干扰能力强，在当前实验 力学领域具有广泛的应用。它的基本原理是对变形 前后的图像子区进行相关性计算，图像子区一般选 取以待求点为中心点的正方形区域。对变形后的子 区按恰当的相关函数进行计算，选取相关性最大的 变形后子区为目标图像子区，来计算图像的全场应 变。选用标准化协方差相关函数[13]，该式为
基金项目：国家自然科学基金(11472145；51008166)
收稿日期：2014-10-30
修回日期：2015-03-11
第一作者简介：陈凡秀，女，1979 年生，博士，青岛理工大学，副教授；研究方向——光测力学的理论与应用。 E-mail：mecfx@
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
246
应用力学学报
第 32 卷
展开如下
α
β
⎡% 0 " 0 "⎤ ⎧ # ⎫ ⎧ #
⎫
p
⎢ ⎢
0
⎢#
Kα 0
0 Kβ %0
0 #
⎥ ⎥ ⎥
⎪ ⎪⎪ ⎨
fα #
⎪ ⎪⎪ ⎬
=
⎪⎪⎪b p ⎨
−K #
β
f
β
⎪ ⎪⎪ ⎬
q⎢⎢" −K α "
0
"⎥⎥
⎪ ⎪
f
β
⎪ ⎪
⎪ ⎪
bq
⎪ ⎪
⎢⎣ # 0 # 0 %⎥⎦ ⎪⎩ # ⎪⎭ ⎪⎩ #
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 为子区中心的位移，以像素为单位。
为接触点到原点的坐标矢量。如图 1 所示，与 P 颗
粒接触的 Ncp 个接触力在颗粒 P 的面积 Ωp 范围内满
足力的平衡，这些接触力包括颗粒与颗粒之间的接
触力 f α 、已知的边界力 f β 和未知的边界力 f γ 。
xα 表示第 α 个接触点到图像坐标原点的坐标矢量， 力矩平衡即 P 颗粒上的所有接触力到图像原点的力 矩满足平衡。
文章编号：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
1引言
在工程实际中，研究颗粒物质的力学性质具 有非常重要的意义。颗粒间力的传播路径（即力 链），决定着颗粒物质宏观及微观的力学性质。近 年来，研究人员通过各种方法对颗粒受力进行研 究。文献[1]研究了二维光弹颗粒体系力链长度的 分布，并得到力链长度随施加力的变化呈指数衰 减的规律。文献[2]使用光弹性测量技术将二维颗粒 中的受力及力的传递过程可视化并对力的大小进行 定量测量。文献[3]用离散元模拟对二维颗粒材料的 受力变形进行了研究，得到内部的力链屈曲引起宏 观上的剪胀现象的结论。文献[4]通过模拟实验研究
⎪⎭
(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
⎪ ⎭
对每个颗粒都可建立力的平衡、力矩平衡、平
移，并获得εx、εy、γxy。对图像进行边缘提取得到圆 心坐标和半径，对每个颗粒面积区域内的应力求和
4 二维钢珠颗粒受集中力试验
取平均，得到各颗粒的平均应力。
4.1 试验系统
二维钢珠颗粒受集中力的试验系统见图2。 CCD相机采用的是德国Basler产品(scA1600-14fm)。 施加压力的装置为岛津牌万能试验机，型号为
P 颗粒上的平均应力为对 Ωp 面积范围内的应
力值求平均，即
图 1 颗粒边界力和颗粒间接触力 Fig.1 Illustration of particle-particle and particle-boundary contacts
考虑线性动量平衡、发散定理和柯西应力张量
的对称性，公式(4)可写为
∑ ( ) σ p =
1
Ncp
sym
f α ⊗ xα
Ωp α =1
(5)
其中： f α 为待求的颗粒间接触力；符号 ⊗ 表示并
列运算，比如 (α ⊗ β )ijkl = αij βkl 。
在二维平面上，式(2)、式(3)、式(5)可用矩阵形 式表示，即
Ncp
∑ K α ⋅ f α =b p
(6)
α =1
AG-IC50kN。该试验机可实时显示钢条底部在钢条 下压过程中的受力并通过相连的计算机绘制成钢条 底部接触力随压入深度的变化图。试验采用直径为 5mm的钢珠，为了防止钢珠表面反光并区分灰度， 在钢珠表面人工形成散斑。在凹形槽中铺入单层的 钢珠颗粒，用透明有机玻璃覆盖并将有机玻璃固定 于凹形槽上作为观测面。将与万能试验机压力探头 连接的钢条(横截面积：1mm×9mm)以恒定速率压 入钢珠颗粒体系中，CCD相机以2幅/秒的频率采集 体系变形图。
网络出版时间：2015-04-20 11:43 网络出版地址：/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
4.3 颗粒体系的排列结构
由于钢珠上喷涂了散斑，直接对原始图片提取 钢珠边缘较难实现，需要对图像进行初步的处理。 首先对初始图像进行二值化处理，如图 3(b)所示， 目的是将钢珠颗粒与背景更好地区分开。采用自适
第2期
陈凡秀，等：基于 DIC 的颗粒间接触力计算及力链分析
247
应 Canny 算子[14]对二值化后的图像进行初步边缘提 取，Canny 算子具有较好的抗噪性能和低误码率， 提取结果见图 3(c)。对初步提取出钢珠边缘的图像 采用 Hough 圆变换进行检测。Hough 圆变换是将平
图 4 钢珠排列结构 Fig.4 Information on the fabric: particle number, position of
⎧⎪ε ⎪
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+μ）γ