基于QGA-SVM的堆石料离散元细观参数标定模型

杨 杰,马春辉,程 琳,冉 蠡,黄志鸿

(1.西安理工大学水利水电学院,陕西 西安 710048;2.西安理工大学省部共建西北旱区生态水利国家重点实验室,陕西 西安 710048)

随着世界范围内水资源开发力度的不断加大,堆石坝已成为众多高坝大库的推荐坝型,目前正朝向300 m级高坝发展。由于受试验原理、试验环境、试验费用等因素的制约,堆石料室内三轴试验难以全面、准确、真实地反应筑坝材料力学特性,是造成堆石坝运行期沉降超出预期的原因之一[1]。近年来,离散元凭借其可从细观尺度探究宏观现象、物理力学关系明确等优势,在土工试验数值模拟方面得到广泛应用,但离散元仍存在细观参数标定困难等问题,需开展进一步研究。如何快速、准确地标定堆石料离散元细观参数,对研究堆石料材料力学性能、确保水利水电工程安全具有重要意义。

颗粒接触模型及细观参数标定是离散元数值模拟土工试验的重要环节,对模拟结果有着直接、明显的影响,也是制约离散元方法发展的瓶颈之一。由于实际岩土体细观结构的复杂性和当前硏究的局限性,尚没有形成一套完善的力学理论用于描述土体细观与宏观参数之间的定量关系,细观参数标定多数采用反复尝试和人为调试的方式,存在较大的盲目性,且效率低、精度差。当前关于离散元细观与宏观参数间互相影响的研究多集中于定性分析,在定量分析方面的研究相对较少。Yang等[2]研究了岩石平行黏结模型中,颗粒数量、细观参数与宏观参数以及单轴抗压强度的关系;Yoon[3]采用中心合成设计对接触黏结模型的细观参数进行标定,并研究了各细观参数与宏观力学特性的相关性;徐小敏等[4]建立了线性接触模型的颗粒法向刚度、颗粒刚度比等细观弹性常数与颗粒材料宏观弹性常数的经验公式;赵彦国等[5]系统研究了平行黏结模型中细观参数对宏观特性的影响,并提出了细观与宏观特性的理论公式;周喻[6]采用BP神经网络建立了岩土体细观与宏观力学参数间的关系;周博[7]拟合了黏性材料内摩擦角、黏聚力的多元非线性公式,定量地描述细观参数和宏观剪切强度参数的联合关系,并给出了切向与法向黏结强度比的建议值;Cheng等[8]采用连续准蒙特卡洛法对土颗粒细观参数进行了标定。在堆石料研究方面,李守巨[9]以堆石料室内三轴试验的实测应力-应变曲线为目标,采用响应面法反演离散元三轴试验模型中的接触刚度、切向刚度、摩擦系数,取得了良好的效果。以上研究为离散元细观参数标定提供了新思路、新方法,但标定模型多对细观参数数目进行简化处理,且模型缺乏处理小样本、多参数、复杂非线性关联的能力;另外上述研究多以岩石的摩尔-库伦强度参数作为宏观参数开展研究,相关成果难以推广至双曲线模型等更为复杂、参数更多的堆石料本构模型中。

近年来,机器学习算法迅速发展,其在处理小样本、非线性、多输出方面具有明显优势,为上述问题的求解提供了新手段。本文针对堆石料离散元三轴试验中存在的细观参数标定困难问题,以室内三轴试验实测应力-应变曲线为标定目标,采用机器学习理论建立细观与宏观参数间关系,以期迅速、准确的完成堆石料细观参数标定。

1 堆石料离散元接触模型

在离散元三轴试验数值模拟过程中,模型的颗粒数目、粒径级配、组装方式、加载速度、接触模型、细观参数等众多因素均会对试件的应力-应变曲线造成影响,其中接触模型及细观参数是影响最大、最难调控的因素。接触模型是离散元描述颗粒变形、受力、运动等状态的基础,离散元接触模型可分为刚度模型、滑动模型、黏结模型三类[6-7]。其中,刚度模型建立了颗粒间接触力和相对位移间的关系;滑动模型建立法向力和切向力之间的关系,判别两个接触体是否发生相对运动;黏结模型是在刚度模型的基础上,明确了法向拉力和切向力的上限值。黏结模型又可分为接触黏结模型和平行黏结模型,当颗粒黏结破坏后将遵循滑动模型。接触黏结模型通过黏结点将球与球连接,能够设定法向与切向黏结力,其原理如图1所示,图中Fn为颗粒的法向接触力,Un为法向位移,kn为法向接触刚度,Fs为颗粒的切向接触力,Fsc为颗粒切向黏结破裂时的接触力,Fsmax为颗粒发生滑动时的接触力,Us为切向位移,ks为切向接触刚度。平行黏结模型假定球体接触后接触处被其他黏结性材料填充,黏结性材料的有效刚度与球体接触点的刚度并联连接,可设定黏结力和黏结力矩。

图1 接触黏结模型原理[7]

根据接触黏结模型基本原理,若颗粒间运动趋势为法向挤压,其压力与位移关系服从线性刚度模型;若颗粒间运动趋势为法向脱离,其法向拉力与法向位移成正比,且法向拉力达到设定法向黏结力时黏结破裂,颗粒不再有法向拉力;若颗粒间运动趋势为切向挤压,其切向力与切向位移成正比,且切向力达到设定切向黏结力时黏结破裂,颗粒将遵循滑动模型。因此,黏结模型是对线性刚度模型的完善与发展。

近年来,国内外关于离散元堆石料三轴试验研究的接触模型选择与细观参数取值统计情况如表1所示,涉及到的细观参数包括:颗粒的法向接触刚度kn,切向接触刚度ks,摩擦系数μ;颗粒的法向黏结力bn,切向黏结力bs;侧墙的法向接触刚度knw1,上下加压板的法向接触刚度knw2,侧墙的切向接触刚度ksw1、上下加压板的切向接触刚度ksw2;孔隙率n等。对于采用接触黏结模型的堆石料离散元三轴试验,通常是颗粒间仍采用线性刚度模型,对于粒径较大的块石由采用接触黏结模型的颗粒簇代替生成,以模拟块石的复杂形状、破碎过程。综合分析可知:①堆石料离散元三轴试验研究所采用的试样,已逐步由简单的二维模型向更为复杂的三维模型发展;②早期的堆石料细观研究多采用简单的线性刚度模型,随着研究的发展,接触黏结模型被广泛用于构建堆石料中的颗粒簇,以深入研究堆石料破碎发展过程等问题;③由于受到众多因素的影响,细观模型参数变化幅度较大甚至存在数量级的差别,难以总结可推广的堆石料细观参数标定准则;④除个别案例外,接触模型若采用接触黏结模型,其接触刚度参数数值较采用线性刚度模型有大幅度减小;⑤虽然接触黏结模型更符合堆石料力学特性,但相比于线性刚度模型,其涉及参数数量更多、取值范围更广、细观与宏观参数间作用机理更为复杂,因此应对接触黏结模型的细观参数确定开展进一步的研究。

表1 堆石料离散元三轴试验接触模型与细观参数取值统计

注：二维试样尺寸为宽度×长度;三维圆柱试样尺寸为直径×高度；三维立方体试样尺寸为宽度×长度×高度。

不同于岩石、土等材料,堆石料由形状不规则、多棱角、排列紧密的颗粒组成,具有咬合力大、抗剪强度高等特点。若仅用圆球模拟堆石料,势必造成颗粒排列形式单一,咬合力较弱。结合上述分析结论,为了更好地模拟堆石料颗粒物理力学特性,提高离散元三轴试验模拟精度,本研究将建立cluster堆石料颗粒簇,颗粒簇将采用接触黏结模型模拟堆石料力学特性。

2 QGA-SVM标定模型建立

2.1 标定模型目标函数

在离散元数值模拟的基础上,建立基于量子遗传算法(quantum genetic algorithm,QGA)和支持向量机(support vector machine, SVM)的堆石料离散元细观参数标定模型,其目标函数为

S(x1,x2,…,xM)=

(1)

式中:x1,x2,…,xM为M个待确定的离散元细观参数;q为应力-应变曲线中提取的定点个数;Fi为QGA-SVM模型的第i个应变值对应的应力值;Ti为室内三轴试验第i个应变值对应的应力实测值。

2.2 QGA-SVM标定模型

QGA是由Narayanan等[21]基于量子基本理论与遗传算法思维建立的一种算法,是目前量子衍生算法中发展最为成熟的算法之一,具有优秀的全局搜索能力。作为机器学习算法发展最热门的算法之一,SVM[22]建立在统计学习理论的Vapnik-chervonenkis(VC)维理论和结构风险最小原理的基础上。SVM可根据有限的样本信息,寻求模型复杂性和学习能力间的最佳折中,具有较强的理论基础,其极值解为全局最优解而非局部最小值,对未知样本有较好的泛化能力。核函数的类型及参数对机器学习算法的性能有至关重要的影响,为进一步提高SVM模型计算能力,本文采用混合核函数[23]:

K(x,xi)=gexp(-‖x-xi‖2/δ2)+

(1-g)(ηxxi+r)d

(2)

式中:K(·)为核函数;x为训练样本的输入;g为组合核函数的待寻优参数;δ为高斯核参数(带宽参数);η、r、d为多项式核参数。

基于QGA-SVM的堆石料离散元细观参数标定模型主要包括训练机器学习模型和搜索细观参数两部分,其计算流程如图2所示。

本文采用拉丁超立方抽样(latin hypercube sampling,LHS)构建离散元细观参数组合,使尽可能少次数的离散元计算能够表达范围更广的细观参数组合。与随机取样法、正交设计法等其他抽样算法相比,LHS具有适用范围广、抽样估值稳定、样本具有更好的代表性和均匀性等优点[24]。

图2 基于QGA-SVM的堆石料细观参数标定模型流程图

凭借强大的表达、并行计算能力,QGA计算性能对参数依赖性小,其取值较为固定,同时SVM核函数的变化范围也较为固定。QGA-SVM细观参数标定模型具有自动适应不同围压、不同试样的优点,尽可能避免传统智能算法需人为修改模型参数、造成分析结果差别较大和模型推广能力差的问题。QGA参数取值和SVM混合核函数参数变化范围如表2所示。

表2 QGA-SVM参数

3 QGA-SVM标定模型验证

为验证模型可行性,采用文献[10,11,25]中的某水电站筑坝石料粒径级配及其室内三轴试验结果进行模型验证,其粒径级配如图3所示。根据粒径级配,模拟试验采用挤压排斥法生成离散元三轴试样,如图4所示。为模拟堆石料应力应变特性,采用接触黏结模型模拟堆石料,其中堆石料颗粒间采用线性刚度模型,对于粒径较大块石由随机颗粒组成的cluster颗粒簇进行代替,cluster颗粒簇内部通过bond键连接。

图3 某水电站筑坝石料粒径级配曲线

图4 离散元三轴试验试样

离散元三轴试件采用圆柱形,尺寸为∅300 mm×650 mm,初始孔隙比为0.35。试件共5 026个颗粒,其中cluster颗粒簇包含508个颗粒。由伺服控制程序控制试验的等压固结、加载、卸载,围压为800 kPa。设定上下加载压盘的运动速度为0.05 m/s,以模拟试样的静力加载。

为实现基于QGA-SVM的堆石料离散元细观参数标定,采用离散元软件PFC(particle flow code)进行离散元三轴试验模拟,使用FISH语言实现细观参数的读入、离散元三轴试验和应力-应变曲线的输出;随后采用Matlab对应力-应变曲线进行插值,以获得0～15%中每0.5%间隔的轴向应变对应应力值,共计31个应力值;最后采用Matlab完成QGA-SVM算法的训练与计算,实现堆石料离散元细观参数标定。

综合分析相关研究,并在一定的试算基础上确定细观参数取值范围,如表3所示。其中侧墙法向刚度约为球体颗粒刚度的1/10,墙体的切向刚度设为0。通过LHS在细观参数取值范围内构建40组参数组合,利用离散元计算相应的应力-应变曲线,并提取其中31个应力值;将40组材料参数组合作为SVM的输入数据,相应应力值作为SVM的输出数据,采用QGA搜索确定SVM核参数,使其性能达到最佳状态;最后,以室内三轴试验值与计算应力值误差最小为目标,发挥SVM快速计算应力-应变曲线的能力,采用QGA全局搜索细观参数。

表3 堆石料离散元三轴试验细观参数取值范围及标定结果

在SVM训练中,QGA计算得SVM性能最佳时的混合核参数组合为:δ=1.377 4、η=0.943 3、r=0.163 6、d=2.227 8、g=0.444 5。SVM与离散元计算值的平均绝对误差(mean absolute error,MAE)为0.015 3,表明SVM对训练数据的拟合精度非常高,训练已完全能够代替离散元实现应力应变计算。随后采用QGA搜索细观参数,最终细观参数标定如表3所示。

通过离散元正算细观参数标定结果,得室内三轴试验、SVM与离散元计算值对比如图5所示。由图5可知:①SVM计算值与室内三轴试验的MAE为0.25,两者在应力-应变曲线的前部与中部吻合度较高,在中前部及后部略有误差,模型整体表现良好,较好地反映了堆石料变形情况;②SVM与离散元正算值的MAE为0.28,两者在应力-应变曲线的前部、后部存在一定误差,表明SVM在前部、后部的训练需加强;③离散元正算值与室内三轴试验的MAE为0.19,两者在应力-应变曲线的前部存在一定误差,在中部、后部拟合效果较好。计算结果中,离散元正算值相较于SVM计算值更为接近室内三轴试验值,表明SVM在拟合应力-应变曲线中仍存在一定误差,致使正算值存在优于SVM计算值的可能。上述分析表明:基于QGA-SVM的堆石料离散元细观参数标定模型是可行的。

图5 各类试验的应变应力对比

为验证QGA优化速度,本文分别采用QGA与遗传算法(GA)标定细观参数,其收敛速度如图6所示。相比于GA,QGA迭代收敛速度快,计算精度高。因此,本文所建立的QGA-SVM标定模型在计算精度、速度方面优势明显。

图6 优化模型收敛速度对比

4 结 语

通过总结国内外关于堆石料细观模型选择与参数选取的研究进展,认为接触黏结模型更适合作为堆石料的离散元细观模型。针对堆石料离散元三轴试验数值模拟中存在的细观参数标定影响因素多、耗时长、成本高等问题,建立基于QGA-SVM的堆石料细观参数标定模型。模型采用QGA算法优化SVM核函数完成机器学习的训练,使其拟合精度、预测精度达到代替离散元计算的要求。随后,模型以SVM计算应力-应变曲线与室内三轴试验实测值差值最小为目标,采用QGA搜索离散元细观参数,完成堆石料离散元三轴试验细观参数标定。通过实例证明,QGA-SVM可快速、精确地标定离散元材料细观模型参数,可将模型进一步扩展至考虑体应变或本构模型参数等方面,具有良好的学术研究、工程应用推广价值。