基于ABAQUS的SBFEM 及在混凝土坝非线性地震响应分析中的应用

陈灯红 刘云龙 林天成 潘子悦

(1.防灾减灾湖北省重点实验室(三峡大学), 湖北 宜昌 443002;2.三峡大学 土木与建筑学院, 湖北 宜昌443002)

我国西南部地区是强震高发区,同时我国大型300 m 级高混凝土坝主要集中在西南部地区,如白鹤滩拱坝的设计地震加速度为0.325g,大岗山拱坝甚至达到了0.557 5g.一旦这些高坝大库发生溃坝灾变,后果不堪设想,因此,地震作用下的高混凝土拱坝抗震安全研究具有重大意义.

在实际工程中,如水工结构抗震领域,大量问题的求解往往伴随着高度的非线性,这些问题被抽象为大规模偏微分方程组,而现有的数学理论往往无法求得解析解,数值算法便应运而生.经典有限元法作为一种成熟的偏微分方程的数值求解方法,已经广泛应用于各类工程问题的求解,各种商用程序中也多采用经典的有限元方法,如商用程序ABAQUS,因其强大的非线性方程组求解能力而被学者和工程师广泛使用.但随着拟建的实际工程体量日益增加、工程服役的环境日益复杂,对工程分析的计算精度、计算效率要求亦相应提高,有必要对现有的工程结构动力分析方法进行进一步的探索.

比例边界有限元法(the scaled boundary finite element,SBFEM)作为近年来的一种新型偏微分方程组求解算法,最早由Wolf、Song[1]提出,以其具有半解析解、降维求解等优势逐渐成为一种主流的数值计算方法.在无限域的理论研究中,比例边界有限元法被广泛用于无限域动力刚度计算[2],这是由于SBFEM 可以将研究域边界向外扩张,进而准确、便捷地表征无限域.在无限域问题的工程应用中,土-结构相互作用[3]、坝体-库水相互作用[4]领域有一大批成果,Qu等[4]提出了基于二维SBFEM 无限水域、半空间无限土体和重力坝结合的数值计算方法.许贺等[5]提出一种耦合SBFEM-FEM 方法用于研究坝体-库水相互作用.Li等[6]将SBFEM 用于研究波在成层无限域的传播问题,并成功应用于计算坝体地基的动力刚度[7].Xu等[8]研究了可压缩库水-坝体的非线性动力相互作用问题.在断裂和损伤问题中,由于SBFEM 的半解析特性,应力强度因子可以直接获得[9],Egger等[10]提出了一种提高应力强度因子计算精度的舒尔解耦和超级收敛恢复技术.Jiang等[11]将SBFEM 和四叉树网格形式结合起来用于研究细观尺度下的混凝土开裂过程.李建波等[12]应用基于水平集算法的扩展比例边界有限元法研究了裂纹扩展问题.刘钧玉等[13]研究了坝体-库水-地基系统的动态断裂问题.Zhang等[14]研究了SBFEM 非局部损伤模型中的应用,在消除网格敏感性的基础上采用四叉树网格加密了局部损伤区域.该方法在三维问题中同样适用[15].杜成斌等[16]研究了SBFEM 在非局部宏观与微观尺度下的损伤模型并进行开裂模拟.Natarajan等[17]将相场模型与SBFEM 结合用于研究脆性材料在动力作用下的断裂.在接触非线性问题中,Zhang和Song[18]提出一种将非匹配网格转化为匹配网格且在局部进行网格细化的技术,增加了将多个独立计算模型结合用于数值仿真的灵活性.Xing等[19]将该方法拓展到三维问题中,采用点对点接触模型对匹配网格的接触面进行离散,降低了求解难度.Zhang等[20]将SBFEM 应用于开裂面的接触分析,并对静动力作用下的裂纹开展问题进行了研究[21].Chen等[22]对结构与半空间外域接触问题进行了研究.

近年来,将SBFEM 与商业软件结合逐渐称为一种主流观点,Ya 等[23]率先对基于ABAQUS 的SBFEM 进行了研究.孔宪京等[24]自主研发了大型岩土工程分析软件,并对高土石坝工程进行了精细化数值分析.但目前与商用程序相结合的研究还有待进一步加强,因此考虑建立一种基于ABAQUS的比例边界有限元工程结构动力分析方法有一定意义.本文在Ya等[23]的研究基础上发展了一种不引入额外自由度的蒙层单元技术,将SBFEM 与ABAQUS中的接触非线性模型结合,建立了考虑横缝接触非线性的坝体-库水-地基相互作用模型.通过与ABAQUS 计算结果对比,说明该方法为工程结构动力分析提供了一种选择.

1 基于ABAQUS的比例边界有限元法

为增强自身的兼容性,商用程序ABAQUS提供了一系列的二次开发接口,其中子程序接口UEL(user defined element)为用户实现自己的理论提供了便捷的方法,文章即是采用该接口将比例边界有限元法嵌入ABAQUS中.

1.1 比例边界有限元法理论基础

比例边界有限元法控制方程为:

系数矩阵为:

式中:[D]为弹性矩阵;ρ为密度;[Jb]为Jacobian矩阵行列式,矩阵[Bi]为形函数矩阵,意义同经典有限元理论.

1.2 比例边界有限元法理论基础用户自定义单元几何拓扑信息

图1展示了有两个单元多面体单元组成的模型图,图2为图1的补充单元拓扑信息.该简单模型的坐标原点为节点3,两个单元均为七面体单元,如果仅仅依据普通的单元节点信息,是无法抽取多面体的所有面节点信息的,补充了额外的拓扑信息以后,通过索引的方式,所有的面节点信息可以很方便地得到.注意单元信息中的符号,该下部单元的第1、5和6个单元的符号为负,以第一个面为例,相对于比例中心坐标(0.5,0.5,0.5),节点1、2、3与比例中心满足右手法则,因此符号为负号,对于其他符号为正的面,它们的面节点编号顺序相对于比例中心不满足右手法则.

图1 多面体单元示例

1.3 HHT时程积分法

HHT(Hilber-Hughes-Taylor)时程积分法是ABAQUS内部采用的一种隐式积分方法,该方法是Newmark积分法的拓展之一.在HHT 时程积分法中多质点体系的运动方程表示为:

1.4 UEL关键量输出规则

ABAQUS与UEL子程序之间传递数据的过程中需要遵循严格的规则.ABAQUS 传进子程序的LFLAGS数组包含了关键的分析类型信息,表明了当前增量步ABAQUS 主程序的分析状态,UEL 计算的关键量即是根据该数组的不同取值而输出到主程序的.UEL 子程序返回给主程序的关键量中最重要的为AMARTIX 矩阵和RHS矩阵,这二者的计算结果根据LFLAGS的值不同而不同.本文主要涉及到LFLAGS(3)和LFLAGS(1)的取值,AMARTIX矩阵和RHS矩阵相应的具体取值如图3所示.

图3 不同分析状态对应的关键矩阵

其中,关键矩阵的不同取值取决于上下标LFLAGS数组的取值,具体的表达式如图4所示.

图4 不同分析状态的关键矩阵取值

1.5 UEL计算流程图

作为成熟的商用程序,ABAQUS有着严格的程序流程设计,采用了 UEL 子程序接口之后,ABAQUS内部的部分计算流程会被替换,自编程序会与ABAQUS主程序有着大量数据交互.图5为基于SBFEM 开发的UEL用户子程序计算流程图,图6为SBFEM_UEL程序中函数的调用关系.

图5 基于UEL的ABAQUS计算流程图

图6 基于SBFEM 的UEL计算流程图

2 耦合SBFEM-FEM 数值实现

文章提出了一种不引入额外自由度的蒙皮单元技术实现了耦合SBFEM-FEM.蒙皮技术是用户自定义单元与ABAQUS丰富的内置功能相结合的一种思路,这种思路的根本目的是利用ABAQUS自带单元库的一些功能以简化用户子程序接口的一些繁琐操作,规避了直接在用户子程序接口中进行大量繁琐的编程,通过直接定义蒙皮单元便可以完成复杂荷载的施加.以下介绍蒙皮技术的两种具体应用方法,分别是全覆盖蒙皮单元技术和通过直接耦合完成用户自定义单元与ABAQUS内置单元库的结合计算.

图7为通过蒙皮单元给用户自定义单元施加重力的示意图,图中等式左端第一项意为用户自定义单元,蓝色为单元节点,为示意绘制了单元的线,实际上ABAQUS是不识别的,左端第二项为ABAQUS 内置单元库中相同几何拓扑信息的单元,等式右端即为蒙皮单元与用户自定义单元结合后的示意图,结合后的单元力学量来自于内置单元和用户自定义单元,颜色为蓝-红相间.

通过蒙皮技术给用户自定义单元施加重力的数学原理见式(4).

式中:下标S表示该力学量来自用户自定义单元,下标F表示该力学量来自经典的有限元法,即ABAQUS的内置单元,下标G 表示重力项,下标ext表示节点的其余迭代不平衡量.

由式(4)可知,方程右端量已经具有了重力项,这是内置单元的形函数插值的结果.这里需注意的是内置单元与用户自定义单元的属性定义的细节,内置单元的属性需要包括密度,这是为了重力的计算与插值,弹性模量与泊松比的数值则为极小值,这就保证内置单元对方程左端的刚度矩阵无贡献,即KF为极小量.而用户自定义单元的属性则只包括弹性模量与泊松比,密度数值则定义为极小量,因为方程左端的质量矩阵项已由MF提供,因此MS为极小值,这里用户自定义单元只提供了刚度矩阵项KS,完成了刚度矩阵项与质量矩阵项的定义后,对于阻尼项,若采用瑞丽阻尼,则由刚度矩阵项与质量矩阵项线性组合而成.

用户自定义单元与ABAQUS 内置单元库结合示意如图8所示.

图8 用户自定单元结合ABAQUS内置单元

用户自定义单元与内置单元的结合在几何上可以简单理解为两种几何体相加的布尔运算.需要注意的是这两种几何体结合的界面不可以存在悬挂节点,否则会导致力学量插值出错.背后的数学原理见式(5)、(6),所有的下标含义同式(4).

上式两种原理表达式是可以互相转化的,不同于给用户自定义单元施加重力,内置单元与外置单元结合的时候二者的计算是互不关联的.同时,二者的属性定义也不同于给用户自定义单元施加重力,这两种单元所有的力学参数都需要正常定义,这种结合方式在丰富了ABAQUS内置单元库的同时即保留了内置单元库参与用户自定义单元计算的可能,又增加了用户自定义单元在ABAQUS中的适用性.

蒙皮技术只是将用户自定义单元与ABAQUS内置功能结合的一种思路,用户自定义单元与内置单元结合参与计算只是其中一个分支,以上只是介绍了其中两种具体实施方法.另外,蒙皮技术的采用也解决了用户自定义单元在ABAQUS中无法可视化的部分问题.用户自定义单元无法在ABAQUS可视化在一定程度上影响了后处理工作,通过解析ODB 文件,将数据格式转化为其他可视化软件可识别的文件是一种思路,但是往往需要大量的编程,而且数据量较大的ODB文件解析非常耗时,因此探索一种简便的用户自定义单元可视化思路是必要的.蒙皮技术通过使用户自定义单元与内置单元共享节点数据在一定程度上达到了位移可视化的目的,在后处理时只需要隐藏用户自定义单元即可.

3 键槽受动荷载啮合接触分析

考虑精细梯形键槽模型,研究键槽受动荷载啮合接触分析,单个键槽平面如图9所示,考虑带有两组键槽的块体几何模型,详细尺寸如图10所示.

图9 单个键槽平面图(单位:cm)

图10 几何模型尺寸图(单位:cm)

如图11所示,边界面S1为加载面,S2为固定约束面,S3为两块体接触面.采用六面体单元对两个块体进行离散,计算模型如图12～13所示,其中图12为协调网格,图13为非协调网格.

图11 加载、边界面示意图

图12 协调网格离散示意图

图13 非协调网格离散示意图

两块体的力学参数见表1.

表1 模型力学参数

在模型边界上施加三角位移荷载,该位移荷载幅值如图14所示,考虑模型尺寸,最终加载面上的位移值为时程图中的2倍,施加的位移量为幅值曲线的2倍,即峰值位移荷载为2 cm,荷载作用时长为10 s,摩擦因数μ=0.6,采用动力隐式算法,增量步时长自适应.

图14 荷载时程(时域)

对该模型进行加载,将位移荷载沿垂直键槽方向施加,即垂直于加载面S1,对加载情况下的模型力学响应进行分析.设置4种工况分别为:

工况1:采用ABAQUS自带的单元,使用协调网格计算;

工况2:采用基于SBFEM 的用户自定义单元,使用协调网格计算,通过蒙皮技术实现接触属性定义;

工况3:采用ABAQUS自带的单元,使用非协调网格计算;

工况4:采用基于SBFEM 的用户自定义单元,使用非协调网格计算,通过蒙皮技术实现接触属性定义.

在接触面上设置3处参考点,如图15所示.

图15 观测点位置示意图

提取参考点的位移时程曲线和接触压力时程曲线如图16所示.发现计算结果与ABAQUS 自带单元吻合较好,采用蒙皮技术后,用户自定义单元的可视化效果与ABAQUS自带单元保持一致.通过该算例验证了耦合SBFEM-FEM 的非线性分析方法可以用于结构动力计算.

图16 观测点动力响应

图17 1 s时刻位移云图

4 耦合SBFEM-FEM 的混凝土坝非线性地震响应分析

4.1 工程概况

以NG5拱坝为研究对象,NG5水电站对南俄河的控制流域达483 km3,占总河流域的2.9%,最大坝高99.00 m,坝顶宽度6.0 m,拱冠梁底宽42.00 m,心线弧长234.838 m,弧高比2.42,拱顶中心角度92.794°,该水电站总库容为3.41 亿m3,死库容为0.63 m3,水库正常蓄水时候的库容量为3.14亿m3,库容系数达34.9%.依据《水工建筑物抗震设计标准》[25],该水利工程按照500年一遇洪水(P=0.2%)涉及,按照10 000年一遇洪水(P=0.01%)校核,水库正常蓄水位高度1 100.00 m.该拱坝平面布置如图18所示.

图18 NG5拱坝平面布置图

4.2 计算模型与材料参数

建立了考虑横缝接触非线性的混凝土坝体-库水-地基动力相互作用模型.其中坝体部分设置一条中缝,几何模型示意如图19所示.采用六面体单元离散,坝体与坝基采用实体C3D8单元模拟,库水采用声学单元AC3D8模拟,计算网格如图20所示,计算模型网格信息见表2.

表2 计算模型网格单元信息

图19 坝体-库水-地基动力相互作用几何模型

图20 NG5拱坝计算模型(单位:m)

模型材料力学参数见表3.采用黏弹性人工边界来模拟远域地基的辐射阻尼效应,以此来考虑坝体-地基相互作用,黏弹性人工边界的物理原件示意如图21所示.对于横缝非线性,采用ABAQUS中的指数接触模型,其模型表达如图22 所示.接触界面采用ABAQUS内部的面-面接触模型定义界面接触属性,接触截面摩擦因数μ=0.6,根据大量试算[26],界面参数c0应取适当微量,p0取值应足够大,那么在动力计算中可保证界面不发生相互贯入现象,文中取c0=5μm,p0=50 GPa.

表3 材料力学参数

图21 黏弹性边界示意图

图22 指数接触模型示意图

4.3 计算荷载

动荷载考虑地震作用,采用人工拟合地震动,峰值加速度为0.93 m/s2,输入地震动的加速度时程曲线如图23所示.

4.4 计算工况

为两种计算模型设计3种工况,见表4,其中工况3为无质量地基对照模型.

表4 计算工况信息

4.5 地震响应结果分析

在地震作用下,对坝体位移、横缝开度以及动水压力等结果进行了对比分析.在坝体关键位置设置若干参考点,参考点的位置如图24所示.

图24 参考点位置

提取了坝体顶部的横缝开度时程如图25所示,3种工况下,横缝开度峰值分别为2.21、2.70、5.27 mm.对比工况1和工况2发现,两组数据吻合良好.并且采用了黏弹性人工边界后,坝体的横缝开度时程响应较无质量地基模型有较大降幅,降幅为58.06%.

图25 点A-B横缝开度时程曲线

绘制坝体坝顶与坝踵部位顺河向相对位移曲线如图26所示,提取顺河向相对位移与相对速度峰值见表5,并给出工况2、3较工况1的增长幅度.

表5 坝顶-坝踵顺河向峰值数据

图26 顺河向相对位移时程曲线

其中,工况3相对位移及相对速度峰值较工况1分别增长了115.44%、141.10%,说明了无质量地基模型相较有质量地基模型会放大结构地震响应.且工况1和工况2的相对位移峰值以及相对速度峰值的相对误差分别为2.40%和2.44%,说明SBFEM 与FEM 两种方法计算结果一致.

考虑坝面左岸和右岸中部位置的两个关键点,提取了关键点部位的动水压力时程曲线,如图27～28所示,并汇总动水压力峰值于表6.

表6 坝面中部关键点动水压力峰值 (单位:Pa)

图27 C点动水压力时程曲线

图28 D 点动水压力时程曲线

对于工况1和工况2,点C、D 二者的动水压力峰值相对误差分别为0.54%、0.42%,结合时程曲线,发现两种工况吻合良好,说明两种算法的计算效果一致.且由表6可知无质量地基模型动水压力峰值较有质量地基模型,具有放大效应.

5 结 论

本文结合比例边界有限元法理论与商用程序ABAQUS的优点,通过用户子程序接口,将SBFEM嵌入到商用程序ABAQUS中.开发了一种不引入额外自由度的蒙皮单元技术.在此基础上以NG5拱坝为研究对象,采用黏弹性动力人工边界条件,考虑接触非线性,建立了耦合比例边界有限元法-经典有限元法的坝体-库水-地基动力相互作用模型,进行非线性地震响应分析,可以得出以下结论:

1)本文提出蒙皮单元技术可以实现比例边界有限元与ABAQUS自带单元、模型等功能相结合,在保留了该理论计算性能的基础上仍能发挥ABAQUS强大的非线性计算性能.

2)采用蒙皮技术后,用户自定义单元的可视化效果与ABAQUS 自带单元保持一致,验证了耦合SBFEM-FEM 的非线性分析方法可以用于结构动力计算.

3)耦合SBFEM-FEM 方法,可以充分考虑大坝-库水-地基的动力相互作用模型、横缝接触非线性以及地震动输入.计算结果与ABAQUS自带单元的结果吻合良好,计算效率高,为复杂工程结构的动力分析提供一种新的方法.