颗粒流方法PFC从细观角度模拟物质结构，实现了对固体介质复杂力学特性与行为特别是岩石破裂物理过程的机制性描述，这是PFC等非连续介质力学方法区别于传统连续介质力学方法的本质所在。

 

不过，非连续介质力学方法需妥善解决因模型中涉及大量接触关系识别与受力运算与分析效率之间的矛盾，为满足工程尺度岩土体变形破坏机制相关深化研究应用的需要，ITASCA非连续-连续介质力学耦合分析技术应运而生，经多年发展形成目前已成熟的FLAC-PFC耦合分析技术。Cai M等（2007）结合日本神奈川水电站地下洞室开挖现场声发射成果与FLAC-PFC耦合分析研究围岩破裂机制或许是该项技术应用的最早案例。

 

本期内容简要介绍如下内容：

1.       ITASCA连续-非连续介质力学耦合分析技术发展历程

2.       早期耦合技术实现原理与特点

3.       新版耦合技术特点与实现原理

4.       应用案例

 

ITASCA连续-非连续介质力学耦合分析技术发展历程

自发布PFC 4.0以来，Itasca即致力于开发形成基于FLAC-PFC平台的耦合计算分析方法。不过，早期方法技术整合度不高（程序耦合程度较低），即分别在FLAC与PFC环境下各自对其负责的模型部分进行求解，然后利用程序的数据接口技术对两个程序模型的耦合部位实现数据的相互传输。早期耦合方法的基本特点是实现方法较为复杂，难以掌握，且某些功能通用程度不高，对使用者的二次开发能力要求也相对较高。自PFC 6.0以来，Itasca实现了FLAC和PFC两个软件的平台整合，即各自可以插件方式运行于彼此的环境中。由此，连续-非连续耦合分析技术的易用性得以极大提升，且耦合功能的通用性也同时得到显著改善。

 

可见，近年来FLAC-PFC耦合分析技术的重点研发方向主要体现在提高软件耦合度和方法通用性两方面，连续-非连续介质力学耦合计算原理则大体沿袭了早期方法，如下图所示。

图1 FLAC-PFC耦合基本原理

 

连续-非连续介质力学耦合计算的基本特点在于：

l   FLAC、PFC各自运行其负责的模型部分（子模型）。其中，FLAC模型部分主要由单元（Zone）构成，而PFC模型部分则通常含有颗粒（Ball）、颗粒团（Clump）等元素；

l   在FLAC、PFC模型运算过程中，连续-非连续介质力学响应同时在模型的耦合边界部位进行数据交换，实现岩石或岩土力学意义上的耦合作用。

l   在计算迭代过程中，速度与力是耦合分析采用的基本耦合变量。其中，由FLAC求解得到的速度通过耦合边界传递至PFC模型部分，由此PFC模型响应得以更新，其对耦合边界部位产生的作用力继而返回至FLAC模型，作为FLAC模型力学响应更新的边界条件。上述迭代过程反复进行，直至FLAC、PFC模型均满足平衡条件，耦合分析终止。

 

耦合边界的定义与识别，耦合变量的解析与传输是连续-非连续介质力学耦合分析技术研发需重点解决的两个环节。

 

早期耦合技术实现原理与特点

PFC 4.0以FISH功能包的方式提供了早期最为基础的连续-非连续介质力学耦合分析功能，由于其实现方法采用FISH明码开发，因此针对功能扩展支持高度灵活的自定义功能，但其明显不足则在于基础功能通用性不强，且其完全基于FISH编辑的操作方式在实际中应用较为复杂，不利于使用者快速掌握。另外，耦合方法的适用范围也仅限于二维分析情形。上述提及的基础功能主要指耦合边界的定义识别和耦合变量的解析与传输。

 

早期基于FISH方法的连续-非连续耦合分析方法可较为自动的用于简单耦合边界的识别与定义，如矩形、圆形等规则形态，如下图所示。

图2 早期耦合计算方法的耦合边界识别的概念

 

上图中，左侧为FLAC程序负责运行的模型部分，其中通过开挖操作形成的矩形孔洞用于容纳PFC模型部分。耦合分析功能包提供了用于耦合边界识别与定义的专用FISH函数，该函数用以识别FLAC、PFC子模型在空间上具有位置叠加（相交）关系的模型元素。据此，将耦合边界定义为FLAC子模型中孔洞边界处单元（Zone）的边（左图中的红色线条，为多段线），以及PFC子模型边界位置处的颗粒（右图中的红色颗粒），以下分别称“耦合线段”和“耦合颗粒”。

 

一旦完成耦合边界的识别与定义，在FLAC、PFC子模型迭代过程中即可通过耦合边界进行耦合变量的数据传输，具体通过程序提供的FISH数据接口功能（Socket）来实现。其中，“耦合线段”和“耦合颗粒”分别接受由对方传递过来的力和速度，力以集中力的方式作用在“耦合线段”的端点部位（实际为FLAC单元Zone的节点），下图以力传输为例说明耦合变量在耦合迭代及数据传输过程中的解译原理。

 图3 早期耦合计算方法的耦合变量传输技术  

 

P表示PFC子模型中某“耦合颗粒”承受的不平衡力，线段0-1表示对应于“耦合线段”的某FLAC单元的边。“耦合线段”的长度为l，“耦合颗粒”与其端点0、1的距离分别为r、r1。相应的，“耦合颗粒”不平衡力P作用在“耦合线段”端点0、1上的集中力分量的矢量依次为F0、F1。由力多边形关系及平衡条件可知：

式中，。“耦合线段”端点速度响应至“耦合颗粒”的解译与传输原理与上述类似，此处不作说明。

 

总体来说，早期耦合技术为基于Itasca软件的连续-非连续耦合分析技术奠定了基础雏形，但在功能细节上存在如下主要不足：

l   耦合边界识别与定义功能仅适用于规则形态，对于复杂模型形态的情形，需采用用户自定义进行必要的功能扩展与干预；

l   对FLAC和PFC而言，可参与耦合分析的模型元素仅分别支持FLAC中的单元（Zone）和PFC中的颗粒（Ball），不支持其他常用模型元素，如FLAC中的结构单元技术，和PFC中颗粒团技术（Clump）等，功能相对单一。

 

以上不足及其软件的整合度不高等架构设计特点均严重制约了早期耦合分析技术的推广和应用。

 

新版耦合技术特点与实现原理

以早期耦合分析方法为基础，若干具体技术环节在新版方法中得到优化升级，在强化问题解决能力与适用性的同时，功能操作实现过程的易用性也得到显著提升。

 

新版技术升级首先表现在对架构设计进行了改造。与早期耦合过程实现时FLAC、PFC需同时运行不同的是，新版Itasca软件将FLAC3D（7.0）和PFC（6.0）的计算和图形引擎设计为通用图形用户界面和循环系统的插件，实现彼此嵌入方式的平台整合，从而为连续-非连续耦合分析技术提供了统一化分析平台；另一重要特点是耦合分析功能的适用性得以显著提升。具体体现在，对耦合边界进一步丰富和具体化，除FLAC中的单元（Zone）和PFC中的颗粒（Ball）外，其他如结构单元、颗粒团等常用模型元素也可参与耦合计算。

1.    耦合技术分类

总体而言，新版耦合技术为适应不同模型元素参与计算设计了如下三种耦合分析方案：

l  一维结构单元耦合（1D Structural Element Coupling）：支持FLAC中的梁（Beam）、桩（Pile）、锚杆/索（Cable）结构单元与PFC中的颗粒（Ball）、颗粒团（Clump）等pebble模型元素进行耦合。耦合变量通过结构单元中的LINK构件进行数据交换传输（Link的概念可详见FLAC3D技术手册），如图4所示：

图4 新版耦合分析技术：一维结构单元耦合

 

l   界面耦合（Interface Coupling或Wall-Zone Coupling）：包括FLAC单元（Zone）表面与二维结构单元（壳Shell、衬砌Liner、土工织物Geogrid等）与PFC中颗粒（Ball）、颗粒团（Clump）等pebble模型元素进行耦合。耦合变量通过在耦合边界部位设置的模型元素Wall来实现数据交换传输，参见图5。图中，耦合边界A、B分别采用了不同的界面耦合技术，前者为FLAC单元Zone表面与PFC耦合，而后者则采用了结构单元与PFC耦合的分析技术。

图5 新版耦合分析技术：界面耦合

l   叠加耦合（Ball-Zone Coupling）：通过FLAC中单元Zone和PFC中颗粒Ball在空间上形成位置叠加关系来实现耦合，如图6所示。

图6 新版耦合分析技术：叠加耦合

 

在以上三种耦合方式中，其中的叠加耦合方法多用于动力分析，旨在通过增加模型元素规模的方式提高动力分析精度，界面耦合分析方法在现实中应用的情形相对更为广泛。

 

2.     耦合分析实现原理

界面耦合方法的基本原理是在可能发生耦合作用的FLAC模型元素的表面（单元Zone的表面或结构单元面）上设置PFC模型组件Wall，并以此作为用于耦合变量交换传递的媒介，实现连续-非连续耦合作用过程，有必要进一步简单说明其理论原理，特别是耦合变量的解译方法。

 

附着设置于FLAC模型元素表面（单元Zone表面或结构单元面）上的PFC模型组件Wall由三角形拼接形成。其中，在Wall生成过程中，若某Zone表面为三角形，则该表面生成一个三角形的Wall，若Zone表面为四边形，则依据该四边形裂化生成两个三角形的WALL。由于FLAC3D中面状结构单元均由三角形构成，因此，附着于其上的WALL是由与结构单元形态一致的三角形构成。

 

图7表示附着于FLAC3D模型元素的一个Wall单元，CP表示位于PFC模型组件上的接触点，同时定义C为位于Wall上对应于CP的接触点，由于接触部位的变形同时存在拉伸、剪切及其扭转作用，因此CP与C两点在空间的位置可能存在不一致的情形。此外，xi表示Wall角点及其附着部位FLAC3D的Zone节点或结构单元节点的坐标，Ai则定义为三个表征三角形的面积（在耦合变量解译过程中，一个WALL单元进一步分解成为3个三角形）。

图7 新版耦合分析技术：耦合变量解译与交换传输

 

同时，定义ri为CP与的xi距离（ri=CP-xi），接触传递至FLAC模型元素的总接触力与总弯矩（PFC中的接触可承受弯矩）分别为F和M，其中M按下式定义：

式中，Mb、（C - CP）×F分别为接触中的弯矩及其因接触位置错位产生的弯矩。

 

在迭代过程中，以上F、M均为已知量。由此，FLAC3D模型元素在WALL附着部位即角点xi处，承受的来自于接触耦合作用CP的集中力Fi和弯矩Mi采用加权法按下式计算得到：

由于ZONE无法直接承受弯矩，当FLAC3D耦合模型元素为Zone时，程序仅解译传输集中力Fi。与此不同，对于FLAC3D耦合模型元素是面状结构单元的情形，集中力Fi和弯矩Mi均能有效传输。

 

应用案例

在丹巴、杨房沟等重点水电工程建设中，ITASCA中国公司综合利用数码摄影测量技术及其FLAC-PFC耦合分析技术建立综合岩体模型，开展各向异性、不连续性、不均匀性及其尺寸效应等复杂岩石力学特性相关拓展性研究。图8为丹巴工程案例，其引水隧洞工程岩体中普遍发育有岩脉，其优势性分布导致岩体呈各向异性、尺寸效应等力学特性。

 

图8 基于FLAC-PFC耦合分析方法的综合岩体模型

 

参考文献：

[1]      王涛, 韩彦辉, 朱永生等. PFC2D/3D颗粒离散元计算方法及应用[M]. 中国建筑工业出版社, 2020.

[2]      PFC 6.0 Manual. Itasca Consulting Group, Inc.

[3]      张春生，刘宁，褚卫江等. 石英云母片岩各向异性特征多尺度精细化描述[J]. 岩石力学与工程学报，2018.