工程科学学报，第40卷，第3期：302-312,2018年3月 Chinese Journal of Engineering,Vol.40,No.3:302-312,March 2018 DOI:10.13374/j.issn2095-9389.2018.03.006:http://journals.ustb.edu.cn 混合颗粒体光弹力链定量提取方法 李飞”四，杨柳”，王金安2，王超” 1)北京科技大学土木与资源工程学院，北京1000832)金属矿山高效开采与安全教有部重点实验室，北京100083 ☒通信作者，E-mail:lifei2016@usth.edu.cn 摘要岩土工程和采矿工程涉及大量的颗粒物质科学和技术难题，定量识别和提取光弹试验颗粒体系的力链网络结构和 分布特征，对于认识和掌握其内部细观力学机理和研究宏观力学行为至关重要.采用彩色梯度均方值(G)算法，建立了不同 粒径的圆形颗粒和方形颗粒的接触力(F)和G的关系；基于数字图像处理技术，提出了识别和区分图像中不同粒径圆形颗粒 和方形颗粒的方法，获得了光弹图片中力链网络结构和力链分布方位.以煤利矿综放开采为实例，对所提出的力链定量提取方 法进行了验证分析，清晰揭示了综放采面矿压形成机理和本质特征.研究表明：单颗粒的F值与G呈单调递增关系，且粒径 越大，F随GC值的增长速度越快：颗粒体系接触力集中分布在0.5F~F(下为平均接触力)，方颗粒、l2mm圆形颗粒平均接触 力较大，多为强力链：中10mm、8mm圆形颗粒的平均接触力较小，多为弱力链.综放开采顶煤和覆岩中力链主要为树状力链， 以竖向发育为主的强力链传递上部主要荷载，横向发育的弱力链对强力链起到侧向支撑的作用.顶煤放出口附近，由于颗粒 侧向移动弱力链丧失，导致强力链减弱甚至消失 关键词光弹试验：混合颗粒；力链；数字图像处理：综放开采 分类号0348 A quantitative extraction method of force chains for composite particles in a photoelastic experiment LI Fei,YANG Liu,WANG Jin-an2,WANG Chao 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)State Key Lab of Education Ministry for High Efficient Mining and Safety of Metal Mines,Beijing 100083,China Corresponding author,E-mail:lifei2016@ustb.edu.en ABSTRACT In geotechnical and mining engineering,numerous particle matters are involved in scientific and technical problems. Recognition and quantification of the structure and distribution of force chain networks using photoelastic experiments,for instance,are significantly important in understanding the internal mechanism of mesoscopic mechanics and studying the macroscopic mechanical be- havior.Based on the algorithm of the mean square value of color gradient ()the correlations of with the contact force (F)of different sizes of round and square particles were established and combined with digital image processing technology.A method was pro- posed for identifying particles in images and distinguishing square particles from round ones,and force chain structures and force chain distributions were obtained in photoelastic images.Using fully mechanized caving mining as an example,the proposed method was veri- fied,and it elaborated the formation and characterization of mining-induced stress in a top-coal caving mining face.The study shows that F monotonically increases with increased G,and larger particle sizes correspond to a faster growth of F.The contact forces of sin- gular particle distribute primarily between 0.5 and 1.5 times average contact force.Also,the average contact forces in square particles and o12mm circular particles are higher and mainly occur in strong force chains,whereas average contact forces in d10mm and d8 mm 收稿日期：2017-0606 基金项目：国家自然科学基金资助项目(U13612034):中国博士后基金资助项目(2017M610047):中央高校基本科研业务费专项基金资助项目 (FRF-TP-16-075A1)
工程科学学报，第 40 卷，第 3 期: 302--312，2018 年 3 月 Chinese Journal of Engineering，Vol． 40，No． 3: 302--312，March 2018 DOI: 10． 13374 /j． issn2095--9389． 2018． 03． 006; http: / /journals． ustb． edu． cn 混合颗粒体光弹力链定量提取方法 李 飞1) ，杨 柳1) ，王金安1，2) ，王 超1) 1) 北京科技大学土木与资源工程学院，北京 100083 2) 金属矿山高效开采与安全教育部重点实验室，北京 100083 通信作者，E-mail: lifei2016@ ustb． edu． cn 摘 要 岩土工程和采矿工程涉及大量的颗粒物质科学和技术难题，定量识别和提取光弹试验颗粒体系的力链网络结构和 分布特征，对于认识和掌握其内部细观力学机理和研究宏观力学行为至关重要． 采用彩色梯度均方值( G2 ) 算法，建立了不同 粒径的圆形颗粒和方形颗粒的接触力( F) 和 G2 的关系; 基于数字图像处理技术，提出了识别和区分图像中不同粒径圆形颗粒 和方形颗粒的方法，获得了光弹图片中力链网络结构和力链分布方位． 以煤矿综放开采为实例，对所提出的力链定量提取方 法进行了验证分析，清晰揭示了综放采面矿压形成机理和本质特征． 研究表明: 单颗粒的 F 值与 G2 呈单调递增关系，且粒径 越大，F 随 G2 值的增长速度越快; 颗粒体系接触力集中分布在 0. 5F ～ F( F 为平均接触力) ，方颗粒、12 mm 圆形颗粒平均接触 力较大，多为强力链; 10 mm、8 mm 圆形颗粒的平均接触力较小，多为弱力链． 综放开采顶煤和覆岩中力链主要为树状力链， 以竖向发育为主的强力链传递上部主要荷载，横向发育的弱力链对强力链起到侧向支撑的作用． 顶煤放出口附近，由于颗粒 侧向移动弱力链丧失，导致强力链减弱甚至消失． 关键词 光弹试验; 混合颗粒; 力链; 数字图像处理; 综放开采 分类号 O348 收稿日期: 2017--06--06 基金项目: 国家自然科学基金资助项目( U13612034) ; 中国博士后基金资助项目( 2017M610047) ; 中央高校基本科研业务费专项基金资助项目 ( FＲF--TP--16--075A1) A quantitative extraction method of force chains for composite particles in a photoelastic experiment LI Fei1) ，YANG Liu1) ，WANG Jin-an1，2) ，WANG Chao1) 1) School of Civil and Ｒesource Engineering，University of Science and Technology Beijing，Beijing 100083，China 2) State Key Lab of Education Ministry for High Efficient Mining and Safety of Metal Mines，Beijing 100083，China Corresponding author，E-mail: lifei2016@ ustb． edu． cn ABSTＲACT In geotechnical and mining engineering，numerous particle matters are involved in scientific and technical problems． Ｒecognition and quantification of the structure and distribution of force chain networks using photoelastic experiments，for instance，are significantly important in understanding the internal mechanism of mesoscopic mechanics and studying the macroscopic mechanical behavior． Based on the algorithm of the mean square value of color gradient ( G2 ) ，the correlations of G2 with the contact force ( F) of different sizes of round and square particles were established and combined with digital image processing technology． A method was proposed for identifying particles in images and distinguishing square particles from round ones，and force chain structures and force chain distributions were obtained in photoelastic images． Using fully mechanized caving mining as an example，the proposed method was verified，and it elaborated the formation and characterization of mining-induced stress in a top-coal caving mining face． The study shows that F monotonically increases with increased G2 ，and larger particle sizes correspond to a faster growth of F． The contact forces of singular particle distribute primarily between 0. 5 and 1. 5 times average contact force． Also，the average contact forces in square particles and 12 mm circular particles are higher and mainly occur in strong force chains，whereas average contact forces in 10 mm and 8 mm
李飞等：混合颗粒体光弹力链定量提取方法 ·303· circular particles are lower and primarily occur in weak force chains.In top-coal caving mining,force chains in top-eoal and overbur- den strata are mainly displayed in tree-like forms.Strong force chains,which extend in a vertical direction,transmit majority overbur- den loads,whereas lateral-developed weak force chains play a role in supporting strong force chains.In the vicinity of the top-coal out- lets,because of the lateral movement of particles toward the mined area,the weak force chains in top coal disappear,resulting in the strong force chains becoming weaker or completely disappearing. KEY WORDS photoelastic experiment:composite particles;force chain:digital image processing:fully mechanized top coal caving mining 在岩土工程中，存在诸多由颗粒物质构成的实 利用偏振光透过具有应力双折射效应的透明材料 际问题，例如流沙、碎屑流，陡岩崩塌问题、尾矿 (如环氧树脂塑料、聚碳酸酯塑料等)，当给材料施 坝边坡破坏问题回以及抗滑桩土拱效应问题间.21 加荷载时，可以获得干涉条纹图、等倾线和等差线 世纪土力学的核心问题是建立土体结构性的数学模 等，通过分析这些条纹信息可以获得材料内部各点 型：从细观角度去研究颗粒体系问题，建立颗粒细观 的应力大小和应力方向.Kramar等、Majmudar 结构与岩土宏观力学性质的联系己经成为当下土力 等回利用光弹试验，研究颗粒物质体系的统计力 学最前沿的研究领域 学，该研究思路是源于分子动力学，将颗粒体系类比 颗粒在自身重力或者外界荷载作用下，会形成 于由分子组成的热力学系统，归纳总结了许多颗粒 强度不一的力链遍布于整个颗粒体系.但实际上只 物质的物理力学特性，获得了令人瞩目的研究成果. 有所占比例较少的强力链作为体系的骨架，起到支 罗格斯大学Kramar等o通过光弹图像提取数学 撑大部分荷载的作用.清华大学孙其诚与金 图论模型，用拓扑学中的同调理论分析力链的结构 峰因提出颗粒一力链一颗粒体系的多尺度结构研 及演化过程，易于分析不同条件下的复杂力链网络 究框架，力链恰好就是连接单个颗粒与颗粒体系的 特征.刘源等采用激光全息光弹法研究受集中 重要桥梁.目前关于力链的研究还停留在唯象研究 荷载作用下的颗粒体系内部力的传递和分布，发现 阶段，主要通过对表征现象进行统计和归纳，得到一 颗粒堆积体内部应力的传递和分布遵循抛物线方 般性规律，但其内部机理尚不明确.例如我们可以 程，不同于经典弹性力学理论的结果.杨荣伟与程 统计力链的分布密度，计算力链强度、方位角等，但 晓辉研制了光弹颗粒直剪仪，发现二元颗粒体系 无法根据边界条件和颗粒材料性质等因素对力链演 的几何结构具有各向异性且平均强度力链出现局部 化趋势做出预判.此外，由于颗粒体系的复杂性和 化特征.郑虎的光弹颗粒直剪试验是基于密度 特殊性，许多关键问题尚未得到解决，比如无法从细 匹配原理，消除了直剪仪底部摩擦力的影响，发现剪 观角度解释在剪切颗粒流中的临界堆积系数小于静 切过程中力链的主方向和剪切室的短对角线一致， 态颗粒体系下的临界堆积系数)，甚至对于强力链 并且颗粒在剪切轴上呈现旋涡分布。王金安等的 的判定问题，目前国内外都未形成统一标准. 对煤矿综放开采顶煤与覆岩力链结构及演化进行了 目前国内外研究颗粒体系的方法众多，主要分 光弹试验研究，获得了综放开采过程中力链构型及 为两类，一类是数值模拟分析方法，其中最有代表性 演化规律.然而，对混合颗粒体力链的定量分析需 的是l979年由Cundall提出的适用于岩土工程领域 要开展进一步的研究.鉴于此，本文通过数字图像 的颗粒离散元法，该方法主要用于研究非连续体介 处理技术，识别混合颗粒体系中的不同形状和不同 质的力学行为.另一类是新型试验方法，例如计算 尺寸颗粒，判定颗粒间接触状态，并用彩色梯度均方 机断层X射线扫描技术(CT技术)，X射线具有强 值(G)法测量颗粒间接触力和接触方向，使光弹试 透射和高分辨率等特性，可以探测颗粒材料的三维 验不仅能运用于圆形颗粒材料，还可以拓展到方形 堆积结构.长江科学院程展林等圆使用CT三轴仪 颗粒材料，拓宽了光弹试验的研究与应用范围. 和计算机图像测量分析系统，研究粗粒土试样在受 1光弹试验 力变形过程中的结构动态变化，并对其位置和位移 进行测量，初步探究了粗粒土颗粒的运动规律.然 1.1试验装置 而CT三轴试验无法给出具体颗粒间接触力的大小 本试验采用的是王金安等研制的颗粒体双 和方向，从而无法建立力链与颗粒之间的力学联系. 轴加载双向流动光弹试验机，该试验机顶部和左侧 本文所采用的是光测弹性力学试验方法，它是 分别有竖向加载电动机和水平加载电动机，电动机
李 飞等: 混合颗粒体光弹力链定量提取方法 circular particles are lower and primarily occur in weak force chains． In top-coal caving mining，force chains in top-coal and overburden strata are mainly displayed in tree-like forms． Strong force chains，which extend in a vertical direction，transmit majority overburden loads，whereas lateral-developed weak force chains play a role in supporting strong force chains． In the vicinity of the top-coal outlets，because of the lateral movement of particles toward the mined area，the weak force chains in top coal disappear，resulting in the strong force chains becoming weaker or completely disappearing． KEY WOＲDS photoelastic experiment; composite particles; force chain; digital image processing; fully mechanized top coal caving mining 在岩土工程中，存在诸多由颗粒物质构成的实 际问题，例如流沙、碎屑流，陡岩崩塌问题［1］、尾矿 坝边坡破坏问题［2］以及抗滑桩土拱效应问题［3］． 21 世纪土力学的核心问题是建立土体结构性的数学模 型; 从细观角度去研究颗粒体系问题，建立颗粒细观 结构与岩土宏观力学性质的联系已经成为当下土力 学最前沿的研究领域． 颗粒在自身重力或者外界荷载作用下，会形成 强度不一的力链遍布于整个颗粒体系． 但实际上只 有所占比例较少的强力链作为体系的骨架，起到支 撑大部分荷载的作用［4--5］． 清华大学 孙 其 诚 与 金 峰［6］提出颗粒—力链—颗粒体系的多尺度结构研 究框架，力链恰好就是连接单个颗粒与颗粒体系的 重要桥梁． 目前关于力链的研究还停留在唯象研究 阶段，主要通过对表征现象进行统计和归纳，得到一 般性规律，但其内部机理尚不明确． 例如我们可以 统计力链的分布密度，计算力链强度、方位角等，但 无法根据边界条件和颗粒材料性质等因素对力链演 化趋势做出预判． 此外，由于颗粒体系的复杂性和 特殊性，许多关键问题尚未得到解决，比如无法从细 观角度解释在剪切颗粒流中的临界堆积系数小于静 态颗粒体系下的临界堆积系数［7］，甚至对于强力链 的判定问题，目前国内外都未形成统一标准． 目前国内外研究颗粒体系的方法众多，主要分 为两类，一类是数值模拟分析方法，其中最有代表性 的是 1979 年由 Cundall 提出的适用于岩土工程领域 的颗粒离散元法，该方法主要用于研究非连续体介 质的力学行为． 另一类是新型试验方法，例如计算 机断层 X 射线扫描技术( CT 技术) ，X 射线具有强 透射和高分辨率等特性，可以探测颗粒材料的三维 堆积结构． 长江科学院程展林等［8］使用 CT 三轴仪 和计算机图像测量分析系统，研究粗粒土试样在受 力变形过程中的结构动态变化，并对其位置和位移 进行测量，初步探究了粗粒土颗粒的运动规律． 然 而 CT 三轴试验无法给出具体颗粒间接触力的大小 和方向，从而无法建立力链与颗粒之间的力学联系． 本文所采用的是光测弹性力学试验方法，它是 利用偏振光透过具有应力双折射效应的透明材料 ( 如环氧树脂塑料、聚碳酸酯塑料等) ，当给材料施 加荷载时，可以获得干涉条纹图、等倾线和等差线 等，通过分析这些条纹信息可以获得材料内部各点 的应 力 大 小 和 应 力 方 向． Kramar 等［5］、Majmudar 等［9］利用光弹试验，研究颗粒物质体系的统计力 学，该研究思路是源于分子动力学，将颗粒体系类比 于由分子组成的热力学系统，归纳总结了许多颗粒 物质的物理力学特性，获得了令人瞩目的研究成果． 罗格斯大学 Kramar 等［10--11］通过光弹图像提取数学 图论模型，用拓扑学中的同调理论分析力链的结构 及演化过程，易于分析不同条件下的复杂力链网络 特征． 刘源等［12］采用激光全息光弹法研究受集中 荷载作用下的颗粒体系内部力的传递和分布，发现 颗粒堆积体内部应力的传递和分布遵循抛物线方 程，不同于经典弹性力学理论的结果． 杨荣伟与程 晓辉［13］研制了光弹颗粒直剪仪，发现二元颗粒体系 的几何结构具有各向异性且平均强度力链出现局部 化特征． 郑虎［14］的光弹颗粒直剪试验是基于密度 匹配原理，消除了直剪仪底部摩擦力的影响，发现剪 切过程中力链的主方向和剪切室的短对角线一致， 并且颗粒在剪切轴上呈现旋涡分布． 王金安等［15］ 对煤矿综放开采顶煤与覆岩力链结构及演化进行了 光弹试验研究，获得了综放开采过程中力链构型及 演化规律． 然而，对混合颗粒体力链的定量分析需 要开展进一步的研究． 鉴于此，本文通过数字图像 处理技术，识别混合颗粒体系中的不同形状和不同 尺寸颗粒，判定颗粒间接触状态，并用彩色梯度均方 值( G2 ) 法测量颗粒间接触力和接触方向，使光弹试 验不仅能运用于圆形颗粒材料，还可以拓展到方形 颗粒材料，拓宽了光弹试验的研究与应用范围． 1 光弹试验 1. 1 试验装置 本试验采用的是王金安等［16］研制的颗粒体双 轴加载双向流动光弹试验机，该试验机顶部和左侧 分别有竖向加载电动机和水平加载电动机，电动机 · 303 ·
·304 工程科学学报，第40卷，第3期 最大能够产生5kN荷载，前进速度为1~100mm· 加载双向流动光弹试验.试验设备及其使用方法详 min-,最大进程为200mm,可以通过位移控制或荷 见文献16-17]. 载控制来确定加载方式.试样容器由钢制边框和钢 使用设备分别获取无偏振镜和有偏振镜时 化玻璃组成，可沿右侧门轴转动打开，通过顶部3 的图像，如图1所示.之后工作主要分为两个部分， mm宽的缝隙将颗粒放入试样容器内部.把试样容 第一部分是用无偏振镜下的图像获取每个颗粒位置 器的底部和右侧的钢制边框加工成连通槽，各由11 和颗粒间接触关系.第二部分是用有偏振镜下的图 个独立钢板挡住，可根据试验具体需要打开挡板，释 像获取每个颗粒的平均接触力，统计力链分布方位 放颗粒.不同的挡板打开或闭合，由排列组合知，该 角及力链强度，从而直观体现颗粒物质体系内部力 装置至少可以模拟1.32×102种不同工况下的双轴 链构造. 图1光弹试验图像.(a)无偏振镜：(b)有偏振镜 Fig.1 Images of photoelastic experiment:(a)image without polarizer:(b)image with polarizer 1.2颗粒材料 表1 Lexan*9030聚碳酸酯板性质参数 本试验使用的颗粒是用Lexan*9030透明聚碳 Table 1 Parameters of properties for Lexan*9030 PC 酸酯板加工而成的，该材料具有良好的透光性和高 性质参数 数值 测试方法 抗冲击性能，易于加工和成型而不会发生破裂.试 密度/(gcm3) 1.20 IS01183 验颗粒尺寸分为4种，分别是8mm×3mm、10 泊松比 0.38 ASTM-D638 mm×3mm、中12mm×3mm的圆柱体颗粒和10mm 抗拉屈服强度/MPa 60 IS0527 ×10mm×3mm的长方体颗粒.试验时用圆形颗粒 抗拉破坏强度/MPa 70 IS0527 模拟相对松散地层，为了避免圆形颗粒因为规律排 弹性模量/MPa 2300 IS0527 列而形成密致结构（即晶体化现象），因此将大、中、 抗弯屈服强度/MPa 90 IS0178 小3种圆形颗粒按照2：9：5的比例均匀混合后再填 挠曲模量/MPa 2300 IS0178 入试样容器中.此外，颗粒在加工过程中会形成残 维卡软化温度/℃ 145 IS0306 余应力，需要通过退火消除残余应力，不同材质和尺 热变形温度/℃ 127 IS075 寸的颗粒的退火过程一般不同.该材料的退火保温 适用于计算少量颗粒受力状况.本文采用的计算颗 温度一般在118~125℃，低于材料的热变形温度. 粒间平均接触力的方法是彩色梯度均方值(G)法， Lexan*9030透明聚碳酸酯板的性质如表1所 该方法不依据颗粒内部条纹级数的具体分布，而是 示 重点考虑单个颗粒内部光强梯度.对于由较多颗粒 1.3接触力计算方法 构成的体系，只需要精确定位每个颗粒位置并计算 通过条纹级数n和条纹延伸方向可以确定主应 每个像素点的光强，即可以处理多颗粒接触力问题. 力差值及其方向，再结合弹性力学解析方法（如侧 在介绍彩色梯度均方值法之前，先介绍灰度梯 向应变法、二向剪应力差法、斜射法和数解法等)， 度均方值法.该方法定义了灰度梯度均方值的概 可以求得该点的主应力大小及其方向 念：灰度图像每一个像素点(i,)均有一个自己的灰 但实际上，由于这种方法需要逐个颗粒进行计 度值1，，取值范围为0~255，亮度由深至浅，0表示 算，工作量巨大且受精度限制，所以这种方法一般只 黑色，255表示白色.每个像素点的灰度值与周围8
工程科学学报，第 40 卷，第 3 期 最大能够产生 5 kN 荷载，前进速度为 1 ～ 100 mm· min － 1，最大进程为 200 mm，可以通过位移控制或荷 载控制来确定加载方式． 试样容器由钢制边框和钢 化玻璃组成，可沿右侧门轴转动打开，通过顶部 3 mm 宽的缝隙将颗粒放入试样容器内部． 把试样容 器的底部和右侧的钢制边框加工成连通槽，各由 11 个独立钢板挡住，可根据试验具体需要打开挡板，释 放颗粒． 不同的挡板打开或闭合，由排列组合知，该 装置至少可以模拟 1. 32 × 1042种不同工况下的双轴 加载双向流动光弹试验． 试验设备及其使用方法详 见文献［16--17］． 使用设备［16］分别获取无偏振镜和有偏振镜时 的图像，如图 1 所示． 之后工作主要分为两个部分， 第一部分是用无偏振镜下的图像获取每个颗粒位置 和颗粒间接触关系． 第二部分是用有偏振镜下的图 像获取每个颗粒的平均接触力，统计力链分布方位 角及力链强度，从而直观体现颗粒物质体系内部力 链构造． 图 1 光弹试验图像． ( a) 无偏振镜; ( b) 有偏振镜 Fig． 1 Images of photoelastic experiment: ( a) image without polarizer; ( b) image with polarizer 1. 2 颗粒材料 本试验使用的颗粒是用 Lexan* 9030 透明聚碳 酸酯板加工而成的，该材料具有良好的透光性和高 抗冲击性能，易于加工和成型而不会发生破裂． 试 验颗粒尺寸分为 4 种，分别是 8 mm × 3 mm、10 mm × 3 mm、12 mm × 3 mm 的圆柱体颗粒和 10 mm × 10 mm × 3 mm 的长方体颗粒． 试验时用圆形颗粒 模拟相对松散地层，为了避免圆形颗粒因为规律排 列而形成密致结构( 即晶体化现象) ，因此将大、中、 小 3 种圆形颗粒按照 2∶ 9∶ 5的比例均匀混合后再填 入试样容器中． 此外，颗粒在加工过程中会形成残 余应力，需要通过退火消除残余应力，不同材质和尺 寸的颗粒的退火过程一般不同． 该材料的退火保温 温度一般在 118 ～ 125 ℃，低于材料的热变形温度． Lexan* 9030 透明聚碳酸酯板的性质如表 1 所 示． 1. 3 接触力计算方法 通过条纹级数 n 和条纹延伸方向可以确定主应 力差值及其方向，再结合弹性力学解析方法( 如侧 向应变法、二向剪应力差法、斜射法和数解法等) ， 可以求得该点的主应力大小及其方向． 但实际上，由于这种方法需要逐个颗粒进行计 算，工作量巨大且受精度限制，所以这种方法一般只 表 1 Lexan* 9030 聚碳酸酯板性质参数 Table 1 Parameters of properties for Lexan* 9030 PC 性质参数 数值 测试方法 密度/( g·cm － 3 ) 1. 20 ISO1183 泊松比 0. 38 ASTM--D638 抗拉屈服强度/MPa 60 ISO527 抗拉破坏强度/MPa 70 ISO527 弹性模量/MPa 2300 ISO527 抗弯屈服强度/MPa 90 ISO178 挠曲模量/MPa 2300 ISO178 维卡软化温度/℃ 145 ISO306 热变形温度/℃ 127 ISO75 适用于计算少量颗粒受力状况． 本文采用的计算颗 粒间平均接触力的方法是彩色梯度均方值( G2 ) 法， 该方法不依据颗粒内部条纹级数的具体分布，而是 重点考虑单个颗粒内部光强梯度． 对于由较多颗粒 构成的体系，只需要精确定位每个颗粒位置并计算 每个像素点的光强，即可以处理多颗粒接触力问题． 在介绍彩色梯度均方值法之前，先介绍灰度梯 度均方值法． 该方法定义了灰度梯度均方值的概 念: 灰度图像每一个像素点( i，j) 均有一个自己的灰 度值 Ii，j ，取值范围为 0 ～ 255，亮度由深至浅，0 表示 黑色，255 表示白色． 每个像素点的灰度值与周围 8 · 403 ·
李飞等：混合颗粒体光弹力链定量提取方法 ·305· 每个分量取值范围均为0~255，亮度由深至 浅，在三维方向上构成RGB色彩空间.对每个分量依 次代入(1)式中，可以得到1VR,12,1TG2,IVB2, 并将三者相加得出彩色图像的1V1,2,即： 1Vl2=1VR2+IVG12+1VB2(4) 之后求单个颗粒彩色G值方法与处理灰度图 像相同，即将单颗粒中包含的所有像素的1V,2 代入(2)式，即可得到每个颗粒的彩色梯度均方值 i-1 G2. 图2G2计算原理示意图 根据以上对单颗粒彩色梯度均方值G的定义， Fig.2 Schematic plot for principle of G 对单个圆颗粒进行逐级加载，如图3(a)所示，拟合 施加外力F与G的关系，如图4(a)所示.从图中可 个像素点的灰度值梯度即为这个点的灰度梯度均方 值1712，可计算如下： 以看出，F随G的增大而增大，并且颗粒粒径越大， F随G2的变化速率越快. =)'+()广 圆形颗粒加载力F与G值之间的对应关系 如下， 22 25 b8mm颗粒： (1) F=172.77(G2)3+1270.6(G2)2+105.21(G2)+ 将单个颗粒中包含的所有像素点的1I,2累 0.8942 (5) 加求平均，这样就得出了所研究的单颗粒的灰度梯 中10mm颗粒： 度均方值G2: F=22224(G2)3-3309.7(G2)2+403.99(G2)- 2.4599 (6) (2) 台台 中12mm颗粒： 式中，N2为颗粒内部不包括边界的像素点的数量 F=-8324.7(G)3+3788.7(G)2+ 基于以上灰度图像分析，那么彩色图像G2可以 164.33(G2)+1.3713 (7) 定义如下：彩色图像每个像素点(i,j》由R(i,j)、G 对单个方颗粒进行逐级加载，如图3(b)所示， (i,》、B(i,》三个值决定，分别对应红(R)、绿(G)、 确定G2与施加外力F的关系，拟合曲线如图4(b) 蓝(B)3个分量，即 所示.当G在0~0.04时，F与G近似呈线性关系， 「R(ij》1 当G2大于0.04时，F随G2变化速率逐渐减缓. I(i,i)= G(i,j》 (3) 10mm×10mm方形颗粒加载力F与G值之间 LB(i,)」 的对应关系： 2N 10N 20N 30N 40N 48N b 2N 10N 20N 30N 40N 48N 图3颗粒逐级加载光弹图.(a）山10mm圆形颗粒：(b)10mm×10mm方形颗粒 Fig.3 Photoelastic images of particles under different loads:(a)10mm round particle:(b)10 mmx 10 mm square particle
李 飞等: 混合颗粒体光弹力链定量提取方法 图 2 G2计算原理示意图 Fig． 2 Schematic plot for principle of G2 个像素点的灰度值梯度即为这个点的灰度梯度均方 值| Δ Ii，j | 2 ，可计算如下: | Δ Ii，j | 2 = [ ( Ii － 1，j － Ii + 1，j ) 2 2 + ( Ii，j － 1 － Ii，j + 1 ) 2 2 ( + Ii － 1，j + 1 － Ii + 1，j － 1 槡 ) 2 2 2 + ( Ii － 1，j － 1 － Ii + 1，j + 1 槡 ) 2 2 ] 2 4 ( 1) 图 3 颗粒逐级加载光弹图． ( a) 10 mm 圆形颗粒; ( b) 10 mm × 10 mm 方形颗粒 Fig． 3 Photoelastic images of particles under different loads: ( a) 10 mm round particle; ( b) 10 mm × 10 mm square particle 将单个颗粒中包含的所有像素点的 | Δ Ii，j | 2 累 加求平均，这样就得出了所研究的单颗粒的灰度梯 度均方值 G2 : G2 = 1 N2 ∑ N i = 1 ∑ N j = 1 | Δ Ii，j | 2 ( 2) 式中，N2 为颗粒内部不包括边界的像素点的数量． 基于以上灰度图像分析，那么彩色图像 G2 可以 定义如下: 彩色图像每个像素点( i，j) 由 Ｒ( i，j) 、G ( i，j) 、B( i，j) 三个值决定，分别对应红( Ｒ) 、绿( G) 、 蓝( B) 3 个分量，即 I( i，j) = Ｒ( i，j) G( i，j) B( i，j ) ( 3) 每个分量取值范围均为 0 ～ 255，亮度由深至 浅，在三维方向上构成 ＲGB 色彩空间． 对每个分量依 次代入( 1) 式中，可以得到| Δ Ｒi，j | 2 ，| Δ Gi，j | 2 ，| Δ Bi，j | 2 ， 并将三者相加得出彩色图像的| Δ Ii，j | 2 ，即: | Δ Ii，j | 2 = | Δ Ｒi，j | 2 + | Δ Gi，j | 2 + | Δ Bi，j | 2 ( 4) 之后求单个颗粒彩色 G2 值方法与处理灰度图 像相同，即将单颗粒中包含的所有像素的 | Δ Ii，j | 2 代入( 2) 式，即可得到每个颗粒的彩色梯度均方值 G2 ． 根据以上对单颗粒彩色梯度均方值 G2 的定义， 对单个圆颗粒进行逐级加载，如图 3( a) 所示，拟合 施加外力 F 与 G2 的关系，如图 4( a) 所示． 从图中可 以看出，F 随 G2 的增大而增大，并且颗粒粒径越大， F 随 G2 的变化速率越快． 圆形颗粒加载力 F 与 G2 值之间的对应关系 如下， 8 mm 颗粒: F = 172. 77 ( G2 ) 3 + 1270. 6 ( G2 ) 2 + 105. 21( G2 ) + 0. 8942 ( 5) 10 mm 颗粒: F = 22224 ( G2 ) 3 － 3309. 7 ( G2 ) 2 + 403. 99( G2 ) － 2. 4599 ( 6) 12 mm 颗粒: F = － 8324. 7 ( G2 ) 3 + 3788. 7 ( G2 ) 2 + 164. 33( G2 ) + 1. 3713 ( 7) 对单个方颗粒进行逐级加载，如图 3( b) 所示， 确定 G2 与施加外力 F 的关系，拟合曲线如图 4( b) 所示． 当 G2 在 0 ～ 0. 04 时，F 与 G2 近似呈线性关系， 当 G2 大于 0. 04 时，F 随 G2 变化速率逐渐减缓． 10 mm × 10 mm 方形颗粒加载力 F 与 G2 值之间 的对应关系: · 503 ·
·306 工程科学学报，第40卷，第3期 50r ( 50r 40 ◆ 30 35 30 20 ◆◆ 。08mm 。φ10mm 20 ▲012mm 15 0.020.040.060.080.100.120.140.16 0.020.040.060.080.100.120.140.16 图4CG2与平均接触力F拟合曲线.(a)圆形颗粒：(b)方形颗粒 Fig.4 Fitted curves of 2 ts average contact force F:(a)round particle:(b)square particle F=31097(G2)3-9620(G2)2+1042.1(G2)+ 程如下： 1.6166 (8) 将图1(a)由彩色图像转换成灰度图像，并提取 基于以上分析可知，当己知光弹试验中每个颗 图像边缘，如图5(a)所示.再选用适当的阈值对图 粒的G时，可以通过代入式(5)~(8)计算出每个 片进行二值化处理，此时，表示边缘的部分呈现白 颗粒间的平均接触力 色，非边缘部分呈现黑色.阈值的选取主要取决于 试验的光线条件，根据多次试验处理经验，此处阈值 2图片信息提取 取值为0.02.在此过程中，图片已经由灰度图像转 2.1无偏振镜下 换成二值图像 从无偏振镜下的图像中可以提取颗粒位置和颗 对储存图片信息的数据矩阵取反，实际处理效 粒接触信息，将这些信息数字化，并储存在数字矩阵 果即黑色转为白色，白色转换为黑色，如图5(b) 所示 中，最后再以图像的形式显示出来，形成颗粒接触网 络，实现接触信息的可视化过程.具体提取算法过 对图5(b)进行“腐蚀”操作“腐蚀”的目的是 a ● 图5无偏振镜图片处理过程.(a)步骤1：(b)步骤2：(c)步骤3：(d)步骤4 Fig.5 Image processing without polarizer:(a)step 1:(b)step 2:(c)step 3;(d)step 4
工程科学学报，第 40 卷，第 3 期 图 4 G2与平均接触力 F 拟合曲线． ( a) 圆形颗粒; ( b) 方形颗粒 Fig． 4 Fitted curves of G2 vs average contact force F: ( a) round particle; ( b) square particle F = 31097 ( G2 ) 3 － 9620 ( G2 ) 2 + 1042. 1( G2 ) + 1. 6166 ( 8) 基于以上分析可知，当已知光弹试验中每个颗 粒的 G2 时，可以通过代入式( 5) ～ ( 8) 计算出每个 颗粒间的平均接触力． 2 图片信息提取 图 5 无偏振镜图片处理过程． ( a) 步骤 1; ( b) 步骤 2; ( c) 步骤 3; ( d) 步骤 4 Fig． 5 Image processing without polarizer: ( a) step 1; ( b) step 2; ( c) step 3; ( d) step 4 2. 1 无偏振镜下 从无偏振镜下的图像中可以提取颗粒位置和颗 粒接触信息，将这些信息数字化，并储存在数字矩阵 中，最后再以图像的形式显示出来，形成颗粒接触网 络，实现接触信息的可视化过程． 具体提取算法过 程如下: 将图 1( a) 由彩色图像转换成灰度图像，并提取 图像边缘，如图 5( a) 所示． 再选用适当的阈值对图 片进行二值化处理，此时，表示边缘的部分呈现白 色，非边缘部分呈现黑色． 阈值的选取主要取决于 试验的光线条件，根据多次试验处理经验，此处阈值 取值为 0. 02． 在此过程中，图片已经由灰度图像转 换成二值图像． 对储存图片信息的数据矩阵取反，实际处理效 果即黑色转为白色，白色转换为黑色，如图 5 ( b) 所示． 对图 5( b) 进行“腐蚀”操作，“腐蚀”的目的是 · 603 ·
李飞等：混合颗粒体光弹力链定量提取方法 ·307· 使颗粒变小并相互分离，使一些未闭合的颗粒区域 的中心坐标，还需要获取正方形的旋转角度，这个角 达到闭合效果.消除图5(b)中的颗粒内部孤立点， 度的确定可以是正方形任意对角线与X轴的夹角， 并进行“开启”操作，目的是消除细的毛刺并断开白 也可以是正方形任意边长与X轴之间的夹角.此 色缝隙连接的区域，还有平滑区域轮廓的作用，结果 处，给出一种适合光弹试验的正方形位置标定算法. 如图5(c)所示. 图7是从图6中截取的一个方颗粒放大图像.由于 进行了上述步骤后，就可以对图5()进行颗粒 颗粒表面污染物或者曝光过程产生的噪点等原因， 识别.所有闭合区域由3部分构成，包括3种尺寸 实际识别出来的正方形颗粒并非是一个完整的正方 的圆颗粒、方颗粒和颗粒间孔隙.首先，计算每一个 形，导致正方形的顶点不易被识别，所以正方形的边 闭合区域的面积(S)和周长(L).面积小于b8mm 和对角线无法准确获取.这里给出一种适合光弹试 圆颗粒面积和大于中l2mm圆颗粒面积的孔隙可以 验方颗粒识别的方法： 首先被剔除 利用上一步获取的方颗粒中心点A(x。,y。)（实 然后，定义系数m为： 际上是白色区域的中心点，当正方形较完整地被识 m=4ms (9) 别出来时，与真实正方形中心误差小于5%)，计算 在标示的区域内部（即白色部分）所有点P(x,y:) 对于标准圆，m取值为1，对于标准正方形，m 到中心点A(xoya)的距离S, 取值为0.785.一般试验过程中，当m∈(0.89~1] S:=√(x:-xo)2+(y:-y0)2 (10) 时，认为该颗粒是圆形颗粒，再根据其面积判别其尺 取距离A最远点为B(x,y),那么AB两点的 寸；当m∈(0.7~0.89]，认为该颗粒是方形颗粒； 连线可以认为是在正方形的对角线上.令AB连线 当m∈0~0.7]时，则可以认为它是形状奇异的颗 与X轴夹角为0，那么正方形旋转角度为0'=0-π/ 粒孔隙.剔除了颗粒孔隙后，如图5()所示. 4(默认当正方形边长与X、Y轴平行时，旋转角度 实际上，由于颗粒或容器表面污染物的存在，一 0=0). 些颗粒没有被识别，还有个别颗粒空隙不可避免的 被识别为颗粒，这些可以通过人工判别的方法进行 剔除.另外，对于边缘不完整的颗粒，可以通过计算 其与边界距离剔除这些颗粒.处理后图片如图6 所示. 图7方颗粒旋转角度示意图 Fig.7 Schematic view of a square particle 根据以上获取的信息：3种圆颗粒的圆心位置 及直径、正方形中心位置、正方形边长和正方形旋转 角度，可以画出识别后的颗粒图8(a) 对于颗粒间接触问题，本文采用如下方法判定 图6识别出的颗粒 Fig.6 Identified particles 方法： 对于圆心位置和正方形中心位置，可以对识别 V属-)'+-)≤分(d+d)（1) 出的颗粒内部各点坐标求平均获得.而圆形的半径 式中，(x1y)、(x2,y2)为两个颗粒的中心坐标，d1、 和正方形的边长，可以通过像素与实际尺寸换算关 d2分别为两个颗粒的直径，对于方颗粒，d取其边 系获得.唯一需要关注的是正方形旋转角度的确 长.入是放大系数，一方面，这是为了避免由于颗粒 定.对于圆形而言，由于圆形的高度对称性，在二维 变形、相机精度、加工因素等造成的误差，另一方面， 平面上只需要圆心坐标和直径大小，即可确定一个 为了保证方颗粒与圆颗粒接触可以被较好识别.在 圆形的具体位置.但对于正方形，不仅需要正方形 本文中，入取值为1.2.根据以上判定方法，可以做
李 飞等: 混合颗粒体光弹力链定量提取方法 使颗粒变小并相互分离，使一些未闭合的颗粒区域 达到闭合效果． 消除图 5( b) 中的颗粒内部孤立点， 并进行“开启”操作，目的是消除细的毛刺并断开白 色缝隙连接的区域，还有平滑区域轮廓的作用，结果 如图 5( c) 所示． 进行了上述步骤后，就可以对图 5( c) 进行颗粒 识别． 所有闭合区域由 3 部分构成，包括 3 种尺寸 的圆颗粒、方颗粒和颗粒间孔隙． 首先，计算每一个 闭合区域的面积( S) 和周长( L) ． 面积小于 8 mm 圆颗粒面积和大于 12 mm 圆颗粒面积的孔隙可以 首先被剔除． 然后，定义系数 m 为: m = 4πS L2 ( 9) 对于标准圆，m 取值为 1，对于标准正方形，m 取值为 0. 785． 一般试验过程中，当 m∈( 0. 89 ～ 1］ 时，认为该颗粒是圆形颗粒，再根据其面积判别其尺 寸; 当 m∈( 0. 7 ～ 0. 89］，认为该颗粒是方形颗粒; 当 m∈［0 ～ 0. 7］时，则可以认为它是形状奇异的颗 粒孔隙． 剔除了颗粒孔隙后，如图 5( d) 所示． 实际上，由于颗粒或容器表面污染物的存在，一 些颗粒没有被识别，还有个别颗粒空隙不可避免的 被识别为颗粒，这些可以通过人工判别的方法进行 剔除． 另外，对于边缘不完整的颗粒，可以通过计算 其与边界距离剔除这些颗粒． 处理后图片如图 6 所示． 图 6 识别出的颗粒 Fig． 6 Identified particles 对于圆心位置和正方形中心位置，可以对识别 出的颗粒内部各点坐标求平均获得． 而圆形的半径 和正方形的边长，可以通过像素与实际尺寸换算关 系获得． 唯一需要关注的是正方形旋转角度的确 定． 对于圆形而言，由于圆形的高度对称性，在二维 平面上只需要圆心坐标和直径大小，即可确定一个 圆形的具体位置． 但对于正方形，不仅需要正方形 的中心坐标，还需要获取正方形的旋转角度，这个角 度的确定可以是正方形任意对角线与 X 轴的夹角， 也可以是正方形任意边长与 X 轴之间的夹角． 此 处，给出一种适合光弹试验的正方形位置标定算法． 图 7 是从图 6 中截取的一个方颗粒放大图像． 由于 颗粒表面污染物或者曝光过程产生的噪点等原因， 实际识别出来的正方形颗粒并非是一个完整的正方 形，导致正方形的顶点不易被识别，所以正方形的边 和对角线无法准确获取． 这里给出一种适合光弹试 验方颗粒识别的方法: 利用上一步获取的方颗粒中心点 A( x0，y0 ) ( 实 际上是白色区域的中心点，当正方形较完整地被识 别出来时，与真实正方形中心误差小于 5% ) ，计算 在标示的区域内部( 即白色部分) 所有点 P( xi，yi ) 到中心点 A( x0，y0 ) 的距离 Si， Si = ( xi － x0 ) 2 + ( yi － y0 槡 ) 2 ( 10) 取距离 A 最远点为 B( x1，y1 ) ，那么 AB 两点的 连线可以认为是在正方形的对角线上． 令 AB 连线 与 X 轴夹角为 θ，那么正方形旋转角度为 θ' = θ － π/ 4( 默认当正方形边长与 X、Y 轴平行时，旋转角度 θ' = 0) ． 图 7 方颗粒旋转角度示意图 Fig． 7 Schematic view of a square particle 根据以上获取的信息: 3 种圆颗粒的圆心位置 及直径、正方形中心位置、正方形边长和正方形旋转 角度，可以画出识别后的颗粒图 8( a) ． 对于颗粒间接触问题，本文采用如下方法判定 方法: ( x1 － x2 ) 2 + ( y1 － y2 槡 ) 2 ≤ λ 2 ( d1 + d2 ) ( 11) 式中，( x1，y1 ) 、( x2，y2 ) 为两个颗粒的中心坐标，d1、 d2分别为两个颗粒的直径，对于方颗粒，d 取其边 长． λ 是放大系数，一方面，这是为了避免由于颗粒 变形、相机精度、加工因素等造成的误差，另一方面， 为了保证方颗粒与圆颗粒接触可以被较好识别． 在 本文中，λ 取值为 1. 2． 根据以上判定方法，可以做 · 703 ·
·308· 工程科学学报，第40卷，第3期 图8无偏振镜下的图片信息提取结果.()颗粒位置分布：(b)接触网络 Fig.8 Extracted figures from images without polarizers:(a)distribution of particles:(b)contact networks 出颗粒接触图，如图8(b) 合均匀的颗粒体系中，大颗粒由于其自身接触面积 2.2有偏振镜下 较大，接触点多，更易于承担荷载，位于强力链中，而 将图1(b)中的光弹彩色图像用平均平方彩色 小颗粒更多起到填充孔隙的作用.当倍数达到3倍 梯度方法处理.一个像素点(i,)会得到一个 后，圆形颗粒的竖坐标显示的频率值接近于0，这就 17112,那么一个1×J像素的彩色图像会得出一 说明圆形颗粒受力为3倍平均接触力时，其几乎等 个I×J的梯度矩阵，再将梯度矩阵转换为图像，如 同于最大受力，再往上的颗粒数目则很少.方形颗 图9所示.图中越亮的地方表示该像素点的 粒的接触力频率分布更接近于中12mm圆形颗粒， 1712越大.从图中可以看出，1V1,2大的地方 方形颗粒的平均接触力更大，也是承担颗粒体系中 主要集中在颗粒接触点、颗粒边缘和位于强力链中 荷载的重要组成部分，且方形颗粒接触力频率分布 的颗粒. 相对集中，几乎不存在大于2倍平均接触力的颗粒 结合图8(a)的颗粒位置信息和图9的像素点 这一规律地发现为后面针对力链强弱的划分提供了 1V12信息，可以利用式(2)求得每个颗粒的G2. 一定的参考. 再由式(5)~(8)可得每个颗粒的接触力F. 05k 一。-◆12mm 。◆10mm 4-◆8mm 0.4 -10mmx10mm 03 0.2 0.1 2 4 56 FIF 图9由彩色数字图像转换成的灰度梯度图 图10接触力频率分布图 Fig.9 Grey gradient image from a color digital image Fig.10 Frequeney distribution of contact force 目前，对于强弱力链的定义和判定方法还存在 通过接触力频率分布曲线（图10）和平均接触 些争议.根据接触力频率分布规律，这里将力链 力（表2），可以得知，4种颗粒的接触力频率分布峰 等级分为强力链、次强力链和弱力链.对于强力链 值点对应的接触力均小于该颗粒的平均接触力F, 定义为：当颗粒的接触力F满足F≥F时，则该颗粒 且方颗粒的平均接触力远大于圆形颗粒的平均接触 位于强力链上.对于次强力链定义为：F>F≥0.5F 力.比较3种尺寸的圆形颗粒分布可以发现，12 时，则该颗粒位于次强力链上.对于弱力链定义为： mm的圆形颗粒分布峰值在0.9倍平均接触力，远 F<0.5F时，则该颗粒位于弱力链上或处于不受外 大于另外两种圆形颗粒，且中12mm的圆形颗粒的 界力作用的状态. 平均接触力也高于另外两种圆形颗粒，这说明在混 在图8(a)的基础上，将强力链和次强力链画
工程科学学报，第 40 卷，第 3 期 图 8 无偏振镜下的图片信息提取结果． ( a) 颗粒位置分布; ( b) 接触网络 Fig． 8 Extracted figures from images without polarizers: ( a) distribution of particles; ( b) contact networks 出颗粒接触图，如图 8( b) ． 2. 2 有偏振镜下 将图 1( b) 中的光弹彩色图像用平均平方彩色 梯度方 法 处 理． 一 个 像 素 点 ( i，j) 会 得 到 一 个 | Δ Ii，j | 2 ，那么一个 I × J 像素的彩色图像会得出一 个 I × J 的梯度矩阵，再将梯度矩阵转换为图像，如 图 9 所 示． 图中越亮的地方表示该像素点的 | Δ Ii，j | 2 越大． 从图中可以看出，| Δ Ii，j | 2 大的地方 主要集中在颗粒接触点、颗粒边缘和位于强力链中 的颗粒． 结合图 8( a) 的颗粒位置信息和图 9 的像素点 | Δ Ii，j | 2 信息，可以利用式( 2) 求得每个颗粒的 G2 ． 再由式( 5) ～ ( 8) 可得每个颗粒的接触力 F． 图 9 由彩色数字图像转换成的灰度梯度图 Fig． 9 Grey gradient image from a color digital image 通过接触力频率分布曲线( 图 10) 和平均接触 力( 表 2) ，可以得知，4 种颗粒的接触力频率分布峰 值点对应的接触力均小于该颗粒的平均接触力 F， 且方颗粒的平均接触力远大于圆形颗粒的平均接触 力． 比较 3 种尺寸的圆形颗粒分布可以发现，12 mm 的圆形颗粒分布峰值在 0. 9 倍平均接触力，远 大于另外两种圆形颗粒，且 12 mm 的圆形颗粒的 平均接触力也高于另外两种圆形颗粒，这说明在混 合均匀的颗粒体系中，大颗粒由于其自身接触面积 较大，接触点多，更易于承担荷载，位于强力链中，而 小颗粒更多起到填充孔隙的作用． 当倍数达到 3 倍 后，圆形颗粒的竖坐标显示的频率值接近于 0，这就 说明圆形颗粒受力为 3 倍平均接触力时，其几乎等 同于最大受力，再往上的颗粒数目则很少． 方形颗 粒的接触力频率分布更接近于 12 mm 圆形颗粒， 方形颗粒的平均接触力更大，也是承担颗粒体系中 荷载的重要组成部分，且方形颗粒接触力频率分布 相对集中，几乎不存在大于 2 倍平均接触力的颗粒． 这一规律地发现为后面针对力链强弱的划分提供了 一定的参考． 图 10 接触力频率分布图 Fig． 10 Frequency distribution of contact force 目前，对于强弱力链的定义和判定方法还存在 一些争议． 根据接触力频率分布规律，这里将力链 等级分为强力链、次强力链和弱力链． 对于强力链 定义为: 当颗粒的接触力 F 满足 F≥F 时，则该颗粒 位于强力链上． 对于次强力链定义为: F ＞ F≥0. 5F 时，则该颗粒位于次强力链上． 对于弱力链定义为: F ＜ 0. 5F 时，则该颗粒位于弱力链上或处于不受外 界力作用的状态． 在图 8( a) 的基础上，将强力链和次强力链画 · 803 ·
李飞等：混合颗粒体光弹力链定量提取方法 ·309· 出，如图11所示.图中红色线表示强力链，蓝色线 表示次强力链.在此基础上，开展统计强力链方位 角以及图8(b)中颗粒接触的方位角，计算强力链、 次强力链所占比例和不同区域颗粒填充分数等一系 列工作，为之后的力链演化分析提供基本数据. 表2平均接触力表 Table 2 Average contact force N 612 mm 中10mm 中8mm 10mm×10mm 9.19 8.38 5.91 29.57 754 B5 B4 B3 B11 BI 1一上部岩层：2一粉砂岩（关键层）：3一砂质泥岩：4一直接顶： 5一煤5层 图12工况模型示意图（单位：mm) Fig.12 Plot of schematic model (unit:mm) 图11力链分布图 Fig.11 Force chain network 释放任何颗粒.从图中可以看出力链主要呈树状分 布，顶部力链比底部强且密，说明上部荷载通过颗粒 3以综放开采为背景的光弹试验 间摩擦和挤压向侧向传递，最后由两侧壁承担部分 3.1工况及模型介绍 上部荷载，这也直观解释了“粮仓效应”的原因.在 颗粒体系在外加荷载作用下，颗粒间相互挤压， 两层方颗粒中，内部应力集中明显，但没有形成明显 力在颗粒体系间传递.但实验证明，力在颗粒体系 力链，而是形成片状区域.经过数字图像处理后，提 中并不是均匀传递的，而是由少量力链构成的网络 取的力链图如图13(b)所示.颗粒接触方位角分 支撑整个颗粒体系的重量和外荷载.某种程度上， 布、强力链方位角分布如图13(c)和图13(d)所示. 颗粒间力和力的分布是被少量外界控制参数明确确 从图13（©）中可以看出，颗粒接触在水平方向明显 定.因此研究力链及力链演化具有重要意义. 多于竖直方向，而其他方向的力链微乎其微.对比 以陈家沟煤矿8512工作面为光弹试验模拟对 图13(c)和图13(d),强力链在竖直方向上所占比 象，本次模拟工作面主采煤层为煤5层，平均埋深 例明显增多，这是由于颗粒体系受到上部荷载作用， 545m,煤层厚度为13.32~54.49m,平均厚度27.8 强力链主要承担竖向荷载 m,煤层平均倾角6.6°.工作面分层厚度8.5~16.6 初始状态下的光弹图像处理结果表明：(1)煤 m,平均开采分层煤厚11.5m.实验中，采用12 层在未开采状态下，岩层内部力链呈树状分布，无明 mm、d10mm和b8mm3种颗粒按2：9：5混合，模拟 显拱形结构：(2)关键层上部岩层的树状力链传递 顶煤及上部覆盖层，采用10mm×10mm方形颗粒， 到关键层，由于关键层强度较高，结构相对稳定，能 每3~5块弱胶结拼接成组合梁模拟关键层及直 够以力键的形式承载上部荷载.而其下面岩层同样 接顶. 用竖向力链支撑着关键层. 试验几何相似比为130、应力相似比为140、密 3.3B5口释放颗粒后的结果分析 度相似比为1.07.建立试验模型如图12所示.该 B3口释放颗粒249个，使B3口上部的第二层 工模型为竖向伺服加载，侧向位移控制，底部从B3 方颗粒落入出口此处，封闭B3口.对B4口和B5口 口向B5口依次释放颗粒，每个口释放颗粒直至第 进行同样的操作步骤，分别释放颗粒64个、102个. 二层方颗粒落于底部为止（直接顶垮落） 此时颗粒状态如图14(a)所示，图14(b)为力链网 3.2初始状态下的结果分析 络图，图中红色箭头表示顶煤放出口位置和顶煤放 图13(a)是初始状态下的光弹图像，此时还未 出方向，黑色箭头表示工作面前进方向，可以看出
李 飞等: 混合颗粒体光弹力链定量提取方法 出，如图 11 所示． 图中红色线表示强力链，蓝色线 表示次强力链． 在此基础上，开展统计强力链方位 角以及图 8( b) 中颗粒接触的方位角，计算强力链、 次强力链所占比例和不同区域颗粒填充分数等一系 列工作，为之后的力链演化分析提供基本数据． 表 2 平均接触力表 Table 2 Average contact force N 12 mm 10 mm 8 mm 10 mm × 10 mm 9. 19 8. 38 5. 91 29. 57 图 11 力链分布图 Fig． 11 Force chain network 3 以综放开采为背景的光弹试验 3. 1 工况及模型介绍 颗粒体系在外加荷载作用下，颗粒间相互挤压， 力在颗粒体系间传递． 但实验证明，力在颗粒体系 中并不是均匀传递的，而是由少量力链构成的网络 支撑整个颗粒体系的重量和外荷载． 某种程度上， 颗粒间力和力的分布是被少量外界控制参数明确确 定． 因此研究力链及力链演化具有重要意义． 以陈家沟煤矿 8512 工作面为光弹试验模拟对 象，本次模拟工作面主采煤层为煤 5 层，平均埋深 545 m，煤层厚度为 13. 32 ～ 54. 49 m，平均厚度 27. 8 m，煤层平均倾角 6. 6°． 工作面分层厚度 8. 5 ～ 16. 6 m，平均开采分层煤厚 11. 5 m． 实验中，采用 12 mm、10 mm 和 8 mm 3 种颗粒按 2∶ 9∶ 5混合，模拟 顶煤及上部覆盖层，采用 10 mm × 10 mm 方形颗粒， 每 3 ～ 5 块弱胶结拼接成组合梁模拟关键层及直 接顶． 试验几何相似比为 130、应力相似比为 140、密 度相似比为 1. 07． 建立试验模型如图 12 所示． 该 工模型为竖向伺服加载，侧向位移控制，底部从 B3 口向 B5 口依次释放颗粒，每个口释放颗粒直至第 二层方颗粒落于底部为止( 直接顶垮落) ． 3. 2 初始状态下的结果分析 图 13( a) 是初始状态下的光弹图像，此时还未 1 ―上部岩层; 2 ―粉砂岩( 关键层) ; 3 ―砂质泥岩; 4 ―直接顶; 5 ―煤 5 层 图 12 工况模型示意图( 单位: mm) Fig． 12 Plot of schematic model( unit: mm) 释放任何颗粒． 从图中可以看出力链主要呈树状分 布，顶部力链比底部强且密，说明上部荷载通过颗粒 间摩擦和挤压向侧向传递，最后由两侧壁承担部分 上部荷载，这也直观解释了“粮仓效应”的原因． 在 两层方颗粒中，内部应力集中明显，但没有形成明显 力链，而是形成片状区域． 经过数字图像处理后，提 取的力链图如图 13 ( b) 所示． 颗粒接触方位角分 布、强力链方位角分布如图 13( c) 和图 13( d) 所示． 从图 13( c) 中可以看出，颗粒接触在水平方向明显 多于竖直方向，而其他方向的力链微乎其微． 对比 图 13( c) 和图 13( d) ，强力链在竖直方向上所占比 例明显增多，这是由于颗粒体系受到上部荷载作用， 强力链主要承担竖向荷载． 初始状态下的光弹图像处理结果表明: ( 1) 煤 层在未开采状态下，岩层内部力链呈树状分布，无明 显拱形结构; ( 2) 关键层上部岩层的树状力链传递 到关键层，由于关键层强度较高，结构相对稳定，能 够以力键的形式承载上部荷载． 而其下面岩层同样 用竖向力链支撑着关键层． 3. 3 B5 口释放颗粒后的结果分析 B3 口释放颗粒 249 个，使 B3 口上部的第二层 方颗粒落入出口此处，封闭 B3 口． 对 B4 口和 B5 口 进行同样的操作步骤，分别释放颗粒 64 个、102 个． 此时颗粒状态如图 14( a) 所示，图 14( b) 为力链网 络图，图中红色箭头表示顶煤放出口位置和顶煤放 出方向，黑色箭头表示工作面前进方向，可以看出， · 903 ·
·310 工程科学学报，第40卷，第3期 (1一上部岩层：2一粉砂岩关键层)： 3一砂质泥岩：4一直接顶：5一煤5层) 90°2000 90°500 120 60 120° 、60 1500 00 150 000 150m 80 180 210 30 210 240 300 240 300 270 270° 图13初始状态.(a)光弹图像：(b)力链网络：(c)颗粒接触方位玫瑰图：(d)强力链方位玫瑰图 Fig.13 Initial state:(a)photoelastic image:(b)force chain networks:(c)rose diagram of particles contact azimuth:(d)rose diagram of strong force chains 在当底部颗粒被释放时，附近的强力链断裂，形成松 压，形成凹形拱结构，即砌体梁结构：(3)由于关键 散结构体系，而在第一层方颗粒上部形成拱状力链， 层的破裂，其上部的岩层力链发生变化，原来的树状 维持上部稳定.较远处的强力链受到扰动，与弱力 力链已经不能直接竖直施加在关键层上，而是以破 链破裂重构，形成新的稳定力链结，同时力链向侧向 裂处为中心，形成凸形拱状力链结构，拱脚位于工作 延伸扩张.原有的方颗粒结构随着底部颗粒的流 面前方和采空区重新压实区域 失，逐渐失稳，以释放口为轴对应倾斜，与水平方向 4结论 形成约15°夹角.颗粒接触方位角分布、强力链方位 角分布如图14(c)和图14(d)所示.对比图13(c) (1)采用彩色梯度均方值的算法，得出不同粒 和图14(c),可以看出颗粒在接触方向上已趋于均 径的圆形颗粒和方形颗粒的F与G值均呈单调递 匀，水平和竖直方向不再占主导地位，整个颗粒体系 增关系，并且颗粒粒径越大，F随G增长速度越快. 向各向同性化演变.再从图14(d)中可知，强力链 (2)用数字图像处理技术，可以有效地消除了 的分布规律基本与颗粒接触分布一致，力链主要集 孔隙和光弹颗粒表面污染物对试验的影响，并且提 中在18°~36°、81°~99°和144°~162°，原有的树状 出识别和区分光弹图像中不同粒径和形状颗粒的算 力链结构己经不足以支撑整个颗粒体系. 法，可以标定颗粒的位置，判断颗粒的接触关系 B5口释放完颗粒后的光弹图像处理结果表明： (3)颗粒体系中，颗粒接触力集中分布在0.5~ (1)放出口周围的力链强度明显减弱甚至消失； 1.5倍平均接触力之间，大的圆形颗粒和方形颗粒 (2)随着煤层开采的推进，煤层上部直接顶逐渐垮 的平均接触力较大，是强力链的重要组成部分，而小 落，直接影响上部关键层结构的稳定性，关键层也逐 的圆形颗粒更多位于次强力链和弱力链上，对强力 渐破裂，但由于岩块排列整齐，岩块转动时会相互挤 链起到侧向支撑的作用
工程科学学报，第 40 卷，第 3 期 图 13 初始状态 . ( a) 光弹图像; ( b) 力链网络; ( c) 颗粒接触方位玫瑰图; ( d) 强力链方位玫瑰图 Fig． 13 Initial state: ( a) photoelastic image; ( b) force chain networks; ( c) rose diagram of particles contact azimuth; ( d) rose diagram of strong force chains 在当底部颗粒被释放时，附近的强力链断裂，形成松 散结构体系，而在第一层方颗粒上部形成拱状力链， 维持上部稳定． 较远处的强力链受到扰动，与弱力 链破裂重构，形成新的稳定力链结，同时力链向侧向 延伸扩张． 原有的方颗粒结构随着底部颗粒的流 失，逐渐失稳，以释放口为轴对应倾斜，与水平方向 形成约 15°夹角． 颗粒接触方位角分布、强力链方位 角分布如图 14( c) 和图 14( d) 所示． 对比图 13( c) 和图 14( c) ，可以看出颗粒在接触方向上已趋于均 匀，水平和竖直方向不再占主导地位，整个颗粒体系 向各向同性化演变． 再从图 14( d) 中可知，强力链 的分布规律基本与颗粒接触分布一致，力链主要集 中在 18° ～ 36°、81° ～ 99°和 144° ～ 162°，原有的树状 力链结构已经不足以支撑整个颗粒体系． B5 口释放完颗粒后的光弹图像处理结果表明: ( 1) 放出口周围的力链强度明显减弱甚至消失; ( 2) 随着煤层开采的推进，煤层上部直接顶逐渐垮 落，直接影响上部关键层结构的稳定性，关键层也逐 渐破裂，但由于岩块排列整齐，岩块转动时会相互挤 压，形成凹形拱结构，即砌体梁结构; ( 3) 由于关键 层的破裂，其上部的岩层力链发生变化，原来的树状 力链已经不能直接竖直施加在关键层上，而是以破 裂处为中心，形成凸形拱状力链结构，拱脚位于工作 面前方和采空区重新压实区域． 4 结论 ( 1) 采用彩色梯度均方值的算法，得出不同粒 径的圆形颗粒和方形颗粒的 F 与 G2 值均呈单调递 增关系，并且颗粒粒径越大，F 随 G2 增长速度越快． ( 2) 用数字图像处理技术，可以有效地消除了 孔隙和光弹颗粒表面污染物对试验的影响，并且提 出识别和区分光弹图像中不同粒径和形状颗粒的算 法，可以标定颗粒的位置，判断颗粒的接触关系． ( 3) 颗粒体系中，颗粒接触力集中分布在 0. 5 ～ 1. 5 倍平均接触力之间，大的圆形颗粒和方形颗粒 的平均接触力较大，是强力链的重要组成部分，而小 的圆形颗粒更多位于次强力链和弱力链上，对强力 链起到侧向支撑的作用． · 013 ·
李飞等：混合颗粒体光弹力链定量提取方法 ·311· ← (1一上部岩层；2一拱结构：3一粉砂岩关键层) 4一砂质泥岩：5一直接顶：6一煤5层：7一采空区) 901000 90300 (e) 1209 09 d 120° 60 800 150 0 150 180 180 210 210 240 300° 240° 300° 270 270 图14B5口释放颗粒后状态.(a)光弹图像：(b)力链网络：(c)颗粒接触方位玫瑰图：(d)强力链方位玫瑰图 Fig.14 State after releasing particles in gate B5.(a)photoelastic image:(b)force chain networks:(c)rose diagram of particles contact azimuth: (d)rose diagram of strong force chains (4)将力链定量提取方法应用于研究综放开采 (张绎，刘洋，吴顺川，等.基于离散一连续耦合的尾矿坝边坡 下顶煤与覆岩的力链结构演化规律，发现覆岩与顶 破坏机理分析.岩土工程学报，2014,36(8)：1473) 煤在初始状态下，力链主要呈竖直方向，上部覆岩荷 B] Han T C,Qiu Z Y,Dou H Q.Soil arching effect between anti- slide piles based on YADE discrete element method.Cent South 载主要通过树状力链形式施加在关键层上，而关键 Unig Sci Technol,2016,47(8):2715 层内部通过力键的形式将荷载传递至顶煤中.随着 (韩同春，邱子义，豆红强.基于颗粒离散元的抗滑桩土拱效 工作面的推进和顶煤的放出，放出口附近的力链减 应分析.中南大学学报（自然科学版），2016,47(8)：2715) 弱甚至消失，关键层上部力链由竖直方向向水平方 [4] Behringer R P,Daniels K E,Majmudar TS,et al.Fluctuations, 向变化，清晰揭示了综放采面矿压形成机理和本质 correlations and transitions in granular materials:statistical me- chanics for a non-conventional system.Phil Trans R Soc London 特征. A,2008,366(1865):493 Kramar M,Goullet A,Kondic L,et al.Persistence of force net- 参考文献 works in compressed granular media.Phys Rev E,2013,87(4): Shi C,Wang S N,Liu L.Research of avalanche disaster numeri- 042207 cal simulation based on granular discrete element method of high- [6]Sun Q C,Jin F.The multi-scale structure of granular matter and steep slope under seismic loads.Chin Rock Mech Eng,2013,32 its mechanics.Physics,2009,38(4):225 (Suppl1):2798 (孙其诚，金蜂.颗粒物质的多尺度结构及其研究框架.物 (石崇，王盛年，刘琳地震作用下陡岩崩塌颗粒离敢元数值 理，2009,38(4)：225) 模拟研究.岩石力学与工程学报，2013,32（增刊1）：2798) ] Wang D M,Zhou Y H.Statistics of contact force network in dense Zhang D.Liu Y,Wu SC,et al.Failure mechanism analysis of granular matter.Particuology,2010,8(2):133 tailing dams based on coupled discrete and continuous method. 8] ChengZL,Wu L P,Ding H S.Research on movement of particle Chin J Geotech Eng,2014,36(8)1473 of fabric of granular material.Rock Soil Mech,2007,28(Suppl):
李 飞等: 混合颗粒体光弹力链定量提取方法 图 14 B5 口释放颗粒后状态 . ( a) 光弹图像; ( b) 力链网络; ( c) 颗粒接触方位玫瑰图; ( d) 强力链方位玫瑰图 Fig． 14 State after releasing particles in gate B5. ( a) photoelastic image; ( b) force chain networks; ( c) rose diagram of particles contact azimuth; ( d) rose diagram of strong force chains ( 4) 将力链定量提取方法应用于研究综放开采 下顶煤与覆岩的力链结构演化规律，发现覆岩与顶 煤在初始状态下，力链主要呈竖直方向，上部覆岩荷 载主要通过树状力链形式施加在关键层上，而关键 层内部通过力键的形式将荷载传递至顶煤中． 随着 工作面的推进和顶煤的放出，放出口附近的力链减 弱甚至消失，关键层上部力链由竖直方向向水平方 向变化，清晰揭示了综放采面矿压形成机理和本质 特征． 参 考 文 献 ［1］ Shi C，Wang S N，Liu L． Ｒesearch of avalanche disaster numerical simulation based on granular discrete element method of highsteep slope under seismic loads． Chin J Ｒock Mech Eng，2013，32 ( Suppl 1) : 2798 ( 石崇，王盛年，刘琳． 地震作用下陡岩崩塌颗粒离散元数值 模拟研究． 岩石力学与工程学报，2013，32( 增刊 1) : 2798) ［2］ Zhang D，Liu Y，Wu S C，et al． Failure mechanism analysis of tailing dams based on coupled discrete and continuous method． Chin J Geotech Eng，2014，36( 8) : 1473 ( 张铎，刘洋，吴顺川，等． 基于离散--连续耦合的尾矿坝边坡 破坏机理分析． 岩土工程学报，2014，36( 8) : 1473) ［3］ Han T C，Qiu Z Y，Dou H Q． Soil arching effect between antislide piles based on YADE discrete element method． J Cent South Univ Sci Technol，2016，47( 8) : 2715 ( 韩同春，邱子义，豆红强． 基于颗粒离散元的抗滑桩土拱效 应分析． 中南大学学报( 自然科学版) ，2016，47( 8) : 2715) ［4］ Behringer Ｒ P，Daniels K E，Majmudar T S，et al． Fluctuations， correlations and transitions in granular materials: statistical mechanics for a non-conventional system． Phil Trans Ｒ Soc London A，2008，366( 1865) : 493 ［5］ Kramar M，Goullet A，Kondic L，et al． Persistence of force networks in compressed granular media． Phys Ｒev E，2013，87( 4) : 042207 ［6］ Sun Q C，Jin F． The multi-scale structure of granular matter and its mechanics． Physics，2009，38( 4) : 225 ( 孙其诚，金峰． 颗粒物质的多尺度结构及其研究框架． 物 理，2009，38( 4) : 225) ［7］ Wang D M，Zhou Y H． Statistics of contact force network in dense granular matter． Particuology，2010，8( 2) : 133 ［8］ Cheng Z L，Wu L P，Ding H S． Ｒesearch on movement of particle of fabric of granular material． Ｒock Soil Mech，2007，28( Suppl) : · 113 ·