一种堆石坝宏观分阶耦合多尺度应力变形计算方法与流程

本发明属于堆石坝数值计算，更具体地，涉及一种堆石坝宏观分阶耦合多尺度应力变形计算方法。

背景技术：

1、据统计，我国已建大坝中土石坝占95％以上，堆石坝作为土石坝的一种，相较于其他坝型，具有施工技术简单、工程造价低、便于就地取材、对地质条件适应性强、抗震性能高等众多优点。坝体沉降、变形是制约堆石坝服役寿命的核心问题之一，超出设计范围的不均匀沉降和变形会致使混凝土面板与堆石体脱空，受力不均加剧，进而产生裂缝，破坏防渗系统，影响堆石坝运行期的安全稳定。

2、近年来，随着现代筑坝设备和技术的发展，堆石坝的设计高度也逐步增加。我国已经有一大批200m、300m高堆石坝陆续投入建设和运行，例如：糯扎渡砾石土心墙堆石坝(261.5m)，雅砻江两河口砾石土心墙堆石坝(295m)，大渡河双江口砾石土心墙堆石坝(312m)，澜沧江如美砾石土心墙堆石坝(315m)。工程实践表明特高堆石坝表现出与中低坝不完全相同的变形性态演变规律，其中一个关键原因在于堆石体天然的多尺度结构特性(细观结构特征会影响堆石体宏观力学特性)在特高堆石坝上体现更为显著。具体而言，由于特高堆石坝的赋存地形地质与运行环境复杂，坝与地基载荷大、应力水平高，以堆石体破碎现象(会造成堆石料的工程设计级配改变，体积收缩变形，进而导致堆石坝产生不均匀沉降和变形)为代表的细观特征更容易发生。

3、堆石坝在细观上由离散的颗粒组成，在宏观上通常处理为连续体。基于有限元等连续介质力学方法，学者们提出了众多的堆石体本构模型，并借助这些本构关系来反映堆石体的非线性、弹塑性、剪胀性和各向异性等，目前在工程实际问题中已得到广泛应用。然而，堆石体应力应变关系以宏观尺度上的经验总结为主，难以从机理上反映堆石体的力学特性。另一方面，细观模型虽能较精确的模拟颗粒组构力学特性过程，把颗粒系统的整个宏观域模型化为单个颗粒的集合体进行数值模拟，但是计算量太大，难以与计算效率达到有机平衡。因此，把细观尺度上的堆石数值模型与宏观尺度上的连续体模型结合起来的多尺度方法则是一种合理的解决方案，能够充分利用连续、非连续介质力学方法这两者的优势，达到计算效率和计算精度的平衡。

4、目前，堆石坝的多尺度计算方法多为层级多尺度法基于均匀化理论，主要包括计算均匀化方法和渐进展开均匀化方法等，一般情况下需要通过表征元来等效宏观单元中积分点的材料性质，最终的问题只在宏观网格上求解，总体来说降低了自由度和计算量，国内外许多学者将该方法包含的计算均匀化方法和渐进展开均匀化方法应用于堆石体的多尺度研究，并取得了丰硕的成果，通常采用有限元法进行宏观分析，但是由于积分点数量庞大，寻找计算精度与效率的平衡仍是亟需解决的问题。一般而言，运行期的高堆石坝的底部与上游区为高应力水平，干湿交替、冻融循环“重灾区”，导致颗粒破碎现象多发，堆石体的细观特征对宏观力学特性影响更为显著，若对该区域采用建立考虑细观结构特性的多尺度模型，其余区域只进行宏观力学分析则可以提升计算效率。

5、本发明的技术基础建立在已授权发明专利“大规模颗粒体破碎模拟分析方法和装置”(公开号：cn112084647a)的基础上，该专利提供了代表体积元的颗粒应力和破碎数值计算理论依据和计算方法。

6、本发明根据宏观有限元模型中的积分点通过颗粒破碎率和应变增量的定性分析，进行分类，对于堆石破碎多发区域，应变增量较高的积分点进行建立多尺度模型，其余区域仍基于有限元法和本构模型计算，建立宏观—分阶耦合多尺度计算方法；对于多尺度分析区域，通过有限元法建立宏观尺度计算模型，利用边界元—离散元联合法进行考虑颗粒破碎和形状的堆石表征元细观尺度计算，并向宏观模型中的积分点反馈应力和切线模量。

技术实现思路

1、针对现有技术的以上缺陷或改进需求，本发明提出了一种堆石坝宏观分阶耦合多尺度应力变形计算方法，具有极高的准确性和有效性。

2、为实现上述目的，按照本发明的一个方面，提供了一种堆石坝宏观分阶耦合多尺度应力变形计算方法，包括：

3、s1：对给定的堆石坝模型进行建模，并利用有限元网格进行数值离散化；

4、s2：结合给定的堆石坝的堆石料的力学参数，建立代表体积元数值模型，进行复杂应力路径下的三轴数值试验，并进行代表体积元的应力、应变与颗粒破碎率的定性分析，以确定宏观有限元模型中的高斯积分点是否需要利用细观代表体积元进行替代的颗粒破碎率阈值；

5、s3：给定一个宏观荷载增量步；

6、s4：求解每个代表体积元的一致性切线模量；

7、s5：结合具体边界条件，利用有限元初步求解整体平衡方程kx＝r，其中，k为总体刚度矩阵；x为节点的待求解位移向量；r为外节点的荷载向量；

8、s6：根据求解得到的每个有限元高斯积分点的变形增量δx，将所有积分点划分为两类：无需代表体积元代替的宏观积分点和需要代表体积元替代计算的细观积分点；

9、s7：针对需要代表体积元替代计算的细观积分点，将变形增量δx作为边界条件进行相应代表体积元的颗粒应力和破碎的模拟计算，直至稳定，将计算得到的应力σ和切线模量d反馈至对应的细观积分点；

10、s8：针对宏观积分点，利用本构模型进行求解得到对应积分点的应力σ和切线模量d；

11、s9：利用公式计算残余力，其中r为残余力，ω为求解域，bt为应变矩阵的转置矩阵，f为外节点荷载向量；

12、s10：重复步骤s4到s9，直至r小于预设的容许残差，则x为最终解；

13、s11：回到步骤s3，继续下一个荷载增量步，直至当前荷载增量步小于最终增量步。

14、在一些可选的实施方案中，步骤s2包括：

15、建立的代表体积元数值模型的边界为周期边界，具体描述为习惯结构沿着如下矢量的平移，其中pα为基矢量，α为正整数，遍历取值1到3。

16、在一些可选的实施方案中，步骤s2包括：

17、代表体积元的应力、应变与颗粒破碎率的定性分析的数值方法为高性能边界元-离散元联合法。所述高性能边界元-离散元联合法的主要步骤包括：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；将颗粒相互作用求解器中某一时间步的计算结果导出，在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件，然后利用边界元法输出待求解颗粒几何模型的内部应力；将每一个颗粒边界上受到的集中力等效为边界上的面力；建立控制方程，并由控制方程得到位移边界积分方程；将位移边界积分方程中体积力导致的域积分转化为边界积分；分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；将每一个颗粒边界上的面力等效为边界单元上的面力；计算单颗粒的系数矩阵，并对系数矩阵求逆；求解边界节点的位移；求解域内的节点的应力和应变等值；根据内部应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中。

18、在一些可选的实施方案中，步骤s2包括：

19、进行代表体积元的应力、应变与颗粒破碎率的定性分析中的应力指的是代表体积元的均匀化应力，利用公式计算得到，其中，σij为平均应力σ的张量形式，np为颗粒个数，v(p)为颗粒p的体积，n为代表体积元内的孔隙率，nc为颗粒p的接触点数量，i和j为遍历取值1到nc的正整数，为颗粒p第i个接触点的坐标位置，为颗粒p形心的坐标位置，颗粒形心指向接触点的单位法向向量分量，为颗粒p第j个接触点的接触力大小。

20、在一些可选的实施方案中，步骤s4包括：

21、求解每个代表体积元的一致性切线模量的方式为：基于均匀应变场的假定，利用以下公式计算，其中v为代表体积元的体积，nc为颗粒间解除数量，kn和ks分别为颗粒间接触点的方向和切向刚度模量，nc和tc分别为颗粒间接触点的外法向和切向单位向量；sc为分支向量。

22、在一些可选的实施方案中，步骤s7包括：

23、位移增量δx需要先被分解为应变增量和旋转增量，再作为边界条件施加。

24、在一些可选的实施方案中，步骤s7包括：

25、进行相应代表体积元的颗粒应力和破碎的模拟计算的数值方法为为高性能边界元-离散元联合法。

26、在一些可选的实施方案中，步骤s8包括：

27、所采用的的本构模型为邓肯-张eb模型。

28、总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

29、根据本发明所提供的一种堆石坝宏观—分阶耦合多尺度应力变形计算方法，相较于传统的宏观有限元模型，能够考虑堆石体细观结构在高应力水平下多发的颗粒破碎现象导致的宏观应力变化，很大程度上可以避免传统有限元模型数值计算过程中出现的“高坝算小，低坝算大”的缺点，提高计算精度。相较于现有的多尺度计算方法，能够将宏观有限元模型中的高斯积分点根据其应变大小划分为宏观和细观积分点，并只针对细观积分点利用代表体积元进行细观应力计算，进而降低了计算量。总体而言，本发明在计算精度和计算效率上取得了很好的平衡，具有可实现性和可靠性，能够应用于堆石坝全尺度的高精度应力分析，以确保工程结构的安全性和可靠性。