基于人工神经网络的颗粒材料本构关系及边值问题研究

张广江, 杨德泽, 楚锡华

(武汉大学 工程力学系, 武汉 430072)

0 引 言

颗粒材料通常是指大量相互接触的离散颗粒与周围流体介质一起组成的复杂体系,广泛存在于日常生产生活的许多领域．以沙石、泥土等典型的颗粒材料为例,其被大量地运用于岩土工程、道路工程、地质工程等工程实践中,还与山体滑坡、土地沙漠化、泥石流等许多地质灾害密切相关．为了保证工程的安全性、预测并防止相关地质灾害的发生,需要提高对颗粒材料的认识,理解并预测其在不同外荷载作用下的力学响应．

目前,研究颗粒材料物理力学性质的主要数值方法有基于离散途径的离散单元法和基于连续体途径的有限单元法．离散单元法[1-2]基于Newton运动定律和力-位移运动定律,通过搜索并计算颗粒之间的相互作用,确定所有颗粒的运动状态．它能够有效地反映颗粒的排列方式、形状、粒度等细观结构信息对颗粒材料力学性质的影响,在计算过程中更加真实地体现颗粒材料的本质,但是由于颗粒数量众多,需要消耗大量的计算资源和时间．有限单元法则是从宏观上将颗粒材料视作一种连续体,采用连续介质理论对其进行分析研究,与离散单元法相比计算效率较高,但分析的结果十分依赖于颗粒材料本构模型的准确性．由于颗粒数量众多、形状各异[3]且内部的细观结构复杂,颗粒材料呈现出非连续性、非均匀性、各向异性[4-5],使其具有了一些十分独特且复杂的物理力学性质,如“粮仓效应”[6-7]、“分离效应”[8-9]、“压力凹陷现象”[10]等,导致建立一个可以精确预测颗粒材料力学行为的本构模型十分困难．现有的颗粒材料本构模型往往是基于唯象假设建立起来的,许多模型十分复杂,本构参数的确定存在一定的困难,或者是仅基于材料的某一现象建立,适用范围存在一定的局限性[11]．

为了解决两者在模拟颗粒材料时的困难,Munjiza和Owen等[12-13]将离散单元法和有限单元法结合,创建了有限元-离散元耦合方法(FDEM),通过有限单元法模拟颗粒材料的变形,内部颗粒的运动则用离散单元法描述,充分发挥了两者的优势．

近年来,随着计算机科学技术的快速发展,算力迅猛提升,新的算法不断涌现,机器学习再次引起了人们的关注．机器学习通过“学习算法”,从现有数据中产生模型,并能通过产生的模型对新的情况进行判断和预测,被认为是一个可能解决材料本构问题的途径[14]．人工神经网络作为机器学习中极具代表性的一类算法,具备强大的学习能力,推动了各相关领域的变革．为此,许多研究人员致力于用其解决颗粒材料数值模拟所面临的一些难题:Benvenuti等[15]利用人工神经网络将宏观实验结果与微观数值参数联系起来,提出了一种识别离散元模拟参数的方法; Zhang 等[16]使用长短期记忆神经网络提出了一种模拟颗粒材料循环特性的方法;Wang等[17]和Qu等[18]分别利用卷积神经网络和循环神经网络成功预测了颗粒材料的本构规律．本文尝试应用人工神经网络算法,将基于离散颗粒模型的离散单元法与基于连续介质模型的有限单元法有机结合,训练一个可以在主空间上准确描述颗粒材料本构关系的模型并运用于求解实际边值问题,形成一套新的、完整的求解方案．

本文基于离散单元法获取颗粒材料的主应力、主应变以及对应的应力-应变矩阵等数据,并利用人工神经网络算法建立起可以在主空间上准确描述颗粒材料本构关系的人工神经网络模型;在假设材料的主应力空间与主应变空间重合的前提下,结合转轴公式,通过UMAT子程序将模型导入ABAQUS中,并用于求解颗粒材料的边值问题．数值算例首先考察了平板双轴压缩问题,结果表明了该方案的可行性和可靠性,并通过进一步考虑边坡问题,结果显示训练后的人工神经网络模型能够在一定程度上预测复杂工况的发展．

1 人工神经网络算法

人工神经网络是由具有适应性的简单单元组成的广泛并行互连的网络,能够模拟生物神经系统对真实世界物体所作出的交互反应[19]．它通常由一个输入层、多个隐藏层和一个输出层组成,且每一层都包含多个神经元,其中输入层通过接收输入值传递给隐藏层,而输出层则产生预测的结果．

从数学上看,人工神经网络本质上是一个复合函数,通过层与层之间的线性变换以及神经元上的非线性激活函数,使其具备了强大的拟合以及预测能力,被广泛地运用于分类以及回归等任务．对于颗粒材料而言,通过唯象假设建立精确的颗粒材料本构模型十分复杂且困难,在已知有关颗粒材料的应力、应变等数据时, 人工神经网络能够直接从数据中学习, 建立应力-应变的内在联系, 十分有效、 便捷地形成颗粒材料的本构模型．

为建立起一个可以准确描述输入数据与输出数据内在关系的人工神经网络模型,需要对人工神经网络不断地进行训练,即通过向前传播和向后传播不断地调整网络中神经元之间的连接权以及阈值,直到损失函数在容许误差范围内．向前传播是指神经网络在训练过程中,数据从输入层传播到输出层的过程．向后传播则是指将神经网络训练后的预测值与训练样本真实值进行比较,将两者之间的误差利用神经网络反向传播调整连接权以及阈值,使得预测值和真实值之间的误差不断下降．

图1 人工神经网络训练过程Fig. 1 The training process of the artificial neural network

2 数据准备与人工神经网络模型的训练

2.1 基于离散元法的数值双轴试验

离散元法是以Newton运动定律为基础、接触力计算为核心,用来求解非连续介质材料问题的一种数值模拟方法,能够有效地对颗粒材料的物理力学行为进行模拟．与基于唯象模型,通过有限元法对颗粒材料的宏观力学行为进行模拟相比,离散元法能够直接把颗粒的排列方式、形状、粒度、孔隙度等细观结构参数以及颗粒的接触本构考虑在内,更加接近颗粒材料的本质属性．为了获取材料在一系列不同外荷载作用下的响应,对于标准试样不论是真实试验还是数值模拟试验,相较于获取其主方向下的应力和应变,获取其一般方向下的应力和应变通常更加困难．因此,本文基于PFC2D,采用离散元法模拟双轴试验,产生颗粒材料的主应力、主应变以及对应的应力-应变矩阵等相关数据,用于训练人工神经网络模型．

离散元数值双轴试验的标准试样如图2所示,其中不同的颜色代表不同的颗粒半径,样本生成模型参数见表1．在一个区域生成约7 200个、半径为2～4 mm随机分布的颗粒后,通过伺服方法,移动试样的四个边界来施加双向等压压力,直到试样稳定达到100 kPa的围压,伺服后通过改变上下、左右墙面的移动速率,来控制双向产生不同的应变比．根据颗粒材料的特性,本文主要考虑以压为主的情况,所有工况的最大轴向应变为7.5%,一共产生了202组数据作为训练集或验证集、40组数据作为测试集,训练集、验证集、测试集的应变加载路径如图3所示．

表1 样本模型参数

图2 离散元数值双轴试样 图3 加载路径示意图 Fig. 2 The discrete element numerical biaxial sample Fig. 3 Schematic diagram of the loading path

为了使训练后的人工神经网络模型能够替代颗粒材料的本构关系用于求解边值问题,离散元数值试验输出的数据主要是主应力、主应变以及对应的应力-应变矩阵．其中,主应力为颗粒样本施加在墙面上的平均反作用力,主应变为墙体的相对位移,应力-应变矩阵[20]为

Cijkl=(kn-ks)aijkl+ksbijkl,

(1)

其中,kn,ks分别为颗粒之间的法向接触刚度以及切向接触刚度;aijkl,bijkl分别为表征元的组构张量．

平面情况下,主应力与主应变的增量关系用矩阵形式表示如下:

(2)

2.2 人工神经网络模型训练过程

为了训练出完整可靠的人工神经网络模型,便于通过ABAQUA子程序UMAT求解实际边值问题,采用了两条神经网络[20]分别训练应变与应力、应力-应变矩阵之间的关系．经过多次试验,人工神经网络结构形式如表2所示．两条神经网络均基于全连接的人工神经网络,采用标准BP算法[21],其中激活函数为sigmoid函数,损失函数为均方误差函数．

表2 神经网络结构参数

为了更快、更精准地训练出人工神经网络模型,所有样本的输入数据、输出数据在进行训练前根据下面公式进行归一化处理,使得所有样本的输入数据、输出数据均位于[-1,1]之间．

(3)

两条神经网络的训练误差如图4所示,为了避免过拟合,在指定训练次数的前提下,取验证集误差最小时的参数为神经网络模型的参数．从图5(a)、(b)可以看出,训练后的人工神经网络模型均能在99%以上解释应变与应力、应力-应变矩阵之间的关系,表明了模型在一定应变范围内具有良好的预测能力．值得注意的是,在后续实际求解边值问题的过程中,训练后的人工神经网络1模型无论是对求解的精度还是求解过程的收敛性都至关重要,而人工神经网络2模型并不影响求解的精度,只对求解过程的收敛性产生重大影响．

(a) 神经网络1 (b) 神经网络2 (a) Neural network 1 (b) Neural network 2

(a) 神经网络1 (b) 神经网络2 (a) Neural network 1 (b) Neural network 2

3 基于人工神经网络模型的UMAT子程序开发

为了将训练后的人工神经网络模型运用于求解边值问题,本文通过UMAT子程序将训练后的人工神经网络模型导入ABAQUS,并用于求解边值问题．以平面问题为例,UMAT子程序流程图如图6所示．

图6 UMAT程序流程Fig. 6 The UMAT program flow chart

为了使得训练后的人工神经网络模型在UMAT中参与计算,假定主应力空间与主应变空间重合,主轴坐标系与一般直角坐标系之间需进行转换．

平面问题转轴公式为:Pij=PmnLimLjn,令σ1应力主方向与x轴之间的夹角为θ．由于主方向上剪应力为0,即

γ12=-2sinθcosθ×εx+2sinθcosθ×εy+(cos2θ-sin2θ)γxy=0,

(4)

得出

(5)

已知θ后,分别通过坐标转换1将应变转换为主应变,坐标转换2将主应力转换为应力:

(6)

在对主方向下的D矩阵进行维度转换后,通过坐标转换3将D矩阵从主应力与主应变的关系转换为一般直角坐标系下应力与应变的关系,即

(7)

(8)

其中

4 数 值 算 例

4.1 双轴压缩试验

为了通过有限元法求解颗粒材料的边值问题,传统方法往往需要将离散元模型中的细观材料参数转换成宏观材料参数,建立以经典弹性模型和Mohr-Coulomb屈服准则为基础形成的弹塑性模型,其函数表达式分别为

σij=Cijklεkl,

(9)

(10)

(11)

其中p为等效应力,q为Mises应力,c为黏聚力,φ为内摩擦角,θ为Lode角．

宏观参数标定方法如下:

1) 在围压达到100 kPa后,保持左右伺服应力不变,在上下两侧缓慢地施加位移荷载,在材料变形为线性变化时,通过下面公式分别计算弹性模量和Poisson比:

(12)

2) 通过改变围压的大小,分别为100 kPa,200 kPa,300 kPa,之后保持左右伺服应力不变,在上下两侧缓慢地施加位移荷载,得到材料的应力-应变曲线(图7(a)),再根据Mohr-Coulomb理论绘制Mohr-Coulomb包络线(图7(b)),并计算相关的材料参数．

(a) 应力-应变曲线 (b) Mohr-Coulomb包络线(a) Stress-strain curves (b) The Mohr-Coulomb envelope

材料的宏观力学参数最终标定为:弹性模量E=6.5×107Pa,Poisson比μ=0.3,内摩擦角32°,剪胀角8°,黏聚力250 Pa．

为了验证离散元法产生的数据训练后的人工神经网络模型是否可以成功替代颗粒材料的本构模型,通过有限元法用于求解实际边值问题,分别运用经典弹塑性模型与人工神经网络模型模拟一个5 m×5 m的平板双轴压缩试验,对比两者的数值计算结果．

离散元模拟过程中通过伺服机制产生了大小为100 kPa的初始应力,使得训练后的人工神经网络模型隐含了该初始应力条件．为了使运用经典弹塑性模型求解与人工神经网络模型求解时的实际加载路径一致,在运用经典弹塑性模型进行求解时,先施加初始应力(图8(a)),再通过位移加载(图8(b))的方式进行双轴压缩,而运用人工神经网络模型求解时,由于训练后的人工神经网络模型隐含了该初始条件,则直接施加相同的位移荷载(图8(b))．

图8 平板压缩试验Fig. 8 The plate compression test

人工神经网络模型与经典弹塑性模型的计算结果如图9所示,分别为位移云图、应变云图和应力云图．在施加相同位移荷载的情况下,应变与应力的计算结果基本一致,其中应变最大值与应力最大值的相对误差分别为9.78%,10.44%．从应变云图可以看出,两者的变形模式相一致,均会出现明显的应变局部化现象,应变主要集中在一个“X”形的剪切带中．为了进一步分析人工神经网络模型与经典弹塑性模型在双轴压缩试验中的计算结果,我们统计了所有单元共计1 600个Gauss积分点处的Mises应力值并计算相对误差,如图10所示,其中柱状图表示位于指定误差范围内的Gauss积分点数目,曲线图表示在指定误差下Gauss积分点数目所占的比例．两个模型计算双轴压缩试验的平均相对误差为11.8%,验证了利用人工神经网络算法建立起可以在主空间上有效描述颗粒材料本构关系的人工神经网络模型,并运用于求解边值问题的可行性．

(a) 经典弹塑性模型(a) The classical elastoplastic model

图10 平板应力误差曲线Fig. 10 The stress error curve of plate

由于经典弹塑性模型和人工神经网络模型的构建方式不同,经典弹塑性模型是基于唯象假设建立的具备一定适用范围的本构模型,在运用于具体材料时,需进行参数标定．而人工神经网络模型则是直接从具体材料的应力、应变等数据中学习,建立了一个适用于该特定材料的本构模型,导致两者在描述材料的本构关系时存在着一定的区别．

经典弹塑性模型存在显式的数学表达式,应力-应变曲线是光滑的,而人工神经网络模型通过数据训练获得,往往不存在显式的数学表达式,其准确性十分依赖于训练数据的数量与质量,而由于离散元数值模拟过程中存在波动导致应力-应变曲线并不是十分光滑,使得训练后的神经网络也不是十分光滑[22]．此外,经典弹塑性模型中弹性模量并不具备应力相关性,即弹性模量始终为一定值,并不会随着材料应力状态而发生改变,而理论上人工神经网络模型能直接从数据中学习,不仅能够反映颗粒细观结构对其宏观力学性能的影响,也能够在一定程度上捕捉材料弹性模量的变化．正是由于经典弹塑性模型和人工神经网络模型本质上存在着这些差异,导致了无论是应力还是应变,两者的计算结果都不可避免地出现一定程度的数值大小差异．

4.2 边坡试验

为了研究训练后的人工神经网络模型在面对更加复杂的情况时,是否具备一定的能力去模拟并预测破坏的发生,进一步基于新建立的求解方案模拟边坡问题．边坡问题是岩土工程中一个十分重要且复杂的问题,不同的边坡在不同的加载条件下会出现不同的破坏形式,其中滑动破坏是最常见的一种破坏类型,为了防止破坏的发生往往需要对边坡的稳定性进行分析．本文模拟边坡问题的数值模型如图11所示,其中边坡顶部为一刚性板,施加偏心位移荷载大小为0.06 m,底部边界为固定约束,右侧边界约束x向位移．图12(a)、(b)分别为经典弹塑性模型和人工神经网络模型的计算结果．

图11 边坡问题示意图Fig. 11 Schematic diagram of the slope

通过对比两者的计算结果,观察到人工神经网络模型与经典弹塑性模型模拟的总体趋势相一致．从应变云图中可以看出,它们均呈现出明显的滑动破坏趋势,且在该位移加载的情况下破坏均未扩展至整个滑坡．通过进一步地提取边坡试验中所有Gauss积分点处的Mises应力并计算相对误差,如图13(a)所示,结果表明,两者计算的相对误差主要集中在50%附近．与平板双轴压缩试验相比,边坡试验由于具备更加复杂的加载情况,其计算误差也相应地明显提高．可见该神经网络模型在面对复杂加载情况时还缺乏足够的能力去精确地捕获材料的应力-应变响应,只能进行定性的分析．

图13 边坡应力误差图Fig. 13 Stress errors curve of slope

与经典弹塑性模型相比,人工神经网络模型的计算误差主要来源于训练误差以及在面对复杂加载情况时,由于主应力空间与主应变空间不重合所造成的误差．虽然训练后的人工神经网络模型在指定的应力-应变范围内能够很好地解释应变-应力、应力-应变矩阵之间的关系,但由于模型泛化能力的限制,图10中仍存在极少数数据点处的相对误差远大于平均误差的情况．而针对主应力空间与主应变空间可能出现不重合所造成的影响,通过提取边坡中应变和Mises应力最大的单元,观察两个模型计算的正应力相对误差随主应力空间与主应变空间夹角变化的趋势,如图13(b)所示,两者总体上呈正相关,正应力的相对误差随着主应力空间与主应变空间夹角的增大而增大,表明了主应力空间与主应变空间的不重合是复杂加载情况下误差的主要来源之一．

5 结 论

本文通过运用人工神经网络算法,将基于离散颗粒模型的离散单元法与基于连续介质模型的有限单元法有机结合,构建了一个可以在主空间上预测颗粒材料本构关系的人工神经网络模型,并运用于求解边值问题,形成了一套完整的求解方案．通过平板双轴压缩试验,验证了该求解方案的可行性,表明了该人工神经网络模型在面对简单加载情况时,能够有效地捕捉应力-应变等数据之间的关系,但是在面对复杂加载情况时,由于主应力空间与主应变空间不重合使得求解结果出现较大的数值差异,导致该方案只能在一定程度上进行定性分析．