基于SPH-FEM耦合方法的落石冲击拱形钢筋混凝土棚洞数值模拟

2019-07-24余志祥郭立平骆丽茹赵世春

振动与冲击 2019年13期

柳春，余志祥,2，郭立平，骆丽茹，赵世春,2

(1. 西南交通大学土木工程学院，成都 610031；2. 陆地交通地质灾害防治技术国家工程实验室，成都 610031)

落石多发于高山峡谷，对铁路和公路沿线造成了巨大潜在威胁[1]。由于落石具有广布偶发和集中群发的特点，主动防护较为困难[2]。棚洞作为一种被广泛应用的落石被动防护结构，主要适用于落石多、不宜采用集中清理的中小型崩塌灾害区域[3]。传统框架钢筋混凝土棚洞一般上覆土垫层，通过土垫层耗散冲击能量(图1(a))，其安全性得到了实际工程的验证，但其结构自重大、跨度小的问题仍然未被解决[4]。目前，利用拱形棚洞能解决上述问题(图1(b))，其结构形式为拱形，不存在大跨度横梁，解决了结构跨度小的问题，另外由于结构形式的改善和轻质缓冲垫层的使用，也解决了自重过大的问题[5]。

目前对钢筋混凝土棚洞动力学行为研究方法有试验方法、理论方法、数值方法。Mougin等[6]等采用试验方法研究了滚石对棚洞混凝土板的冲击。Kishi等[7]对预应力钢筋混凝土棚洞进行了足尺落石冲击测试，以评估棚洞的极限承载能力。Delhomm等[8]通过在混凝土棚洞支座处增设耗能减震器(SDR)以替代砂石垫层吸收落石冲击能量，并进行了缩尺试验研究其动力学行为。Kishi等[9]为了优化冲击力的荷载分配，在钢筋混凝土棚洞上覆三层缓冲层，并进行了足尺试验，最大冲击能量达到了1 500 kJ。Bhatti等[10]进行了拱形棚洞在不同冲击能量下的足尺试验研究。王爽等通过室内模型试验研究了拱形棚洞的抗落石冲击的力学性能，并进行了不同冲击能量、不同冲击位置、不同缓冲层厚度的一些参数分析。但有限的试验并不能代表复杂多变的现场落石冲击条件，若要通过试验全面研究拱棚洞动力学性能，则会造成试验费用较高、周期长。理论研究多集中于落石冲击力研究。刘茂[11]基于弹塑性赫兹接触理论提出了落石冲击力计算公式。易伟等[12]基于量纲分析法提出了改进的落石冲击力计算公式。王林峰等[13]根据消能棚洞的结构特性和受力条件,建立消能棚洞的振动模型,推导消能棚洞刚度系数的计算公式,进而得到消能棚洞的落石冲击力计算公式。国外常用冲击力计算公式基于半经验半理论算法。Montani等[14]以缓冲层的弹性模量、落石速度、质量等为主要考虑因素计算冲击力，但没有考虑冲击体与被冲击体的材料特性，其计算局限性大。日本道路协会根据Kawahara等[15]的研究成果制定了相关落石冲击力规范，但其计算公式的物理量量纲不和谐，其计算结果有一定的偏差。目前，通过理论方法获得物理参量的解析解都存在一定的偏差或地域性。因此数值模拟方法以其经济性与高效性逐渐成为此类冲击问题重要的研究手段，目前主流算法是有限元方法[16-17]，但由于拱棚洞的砂土垫层是由大量离散颗粒体组成的集合体，当受到冲击时，会出现类似流体的大变形性质，若用有限元方法计算，很难处理网格大变形而导致计算失败的问题[18]。而离散元方法存在计算量大和参数标定困难的问题[19]。拱棚洞受冲击过程中存在材料大变形、高速度、高应变率等问题，为研究其动力学行为，则需要一个合理的数值计算方法。

光滑粒子流体动力学方法是一种无网格的连续介质力学计算方法，具有很好的自适应特点，可以很好地处理大变形和后失稳问题，避免了FEM中网格大变形导致计算失败，也避免了DEM中计算量大和参数标定困难等问题。该方法能够描述砂土材料的大变形特征，并具有较高的计算精度和稳定性，但计算时间会有所增加[20-21]。目前，对于钢筋和混凝土的模拟，其材料模型及其有限元算法已比较成熟，使用有限元模拟已能保证足够的精度[22-23]。

(a) 框架棚洞(b) 拱棚洞

图1 钢筋混凝土棚洞

Fig.1 Reinforced concrete tunnel

本文提出了有限元法与光滑粒子流体动力学相耦合的方法研究钢筋混凝土拱棚洞的动力学行为。为提高计算效率和精度，利用SPH粒子模拟落石冲击区域的砂土，非冲击区域砂土用有限元单元模拟；将混凝土、钢筋、岩石、冲击锤等划分为有限元网格。然后将SPH与FEM进行耦合计算。将该方法得到的模拟结果与试验结果进行对比分析，进而为拱形棚洞的研究设计和应用提供一个合理的计算方法。

1 计算理论

1.1 SPH算法基础

SPH方法是一种无网格粒子方法，通过使用一系列任意分布的粒子来求解具有各种边界条件的积分方程或偏微分方程，得到精确稳定的数值解。SPH方程的构造按两个关键步骤进行，第一步为核函数插值，实现场变量或场变量梯度的插值；第二步为粒子近似，实现对核函数估计积分表达式的粒子离散[24]。

在SPH方法中,任意宏观变量f(r)在空间某一点r上的核估计可以通过函数f(r)在定义域中的积分获得

(1)

式中：W(r-r′,h)为核函数，满足正则化、紧支性、非负性、对称性等条件；h为光滑长度，本文取为1.2d0，d0为相邻粒子间距离。

为了能够最终获得粒子的离散控制方程，通过对支持域内一系列粒子的离散化求和来实现，即粒子i处的场函数值可以通过核函数对该粒子紧支域内所有粒子的函数值加权平均得到

(2)

式中：ρj、mj分别为粒子j的密度和质量。

通过式(2)，质量、动量、能量守恒的离散化方程可相应得到。

SPH方法通过使用光滑核函数进行积分表示，它不仅决定了函数近似式的形式、定义粒子支持域的尺寸，还决定了核近似和粒子近似的一致性和精度。本文采用三次样条型核函数来处理三维问题

(3)

1.2 SPH-FEM耦合

在交界面处，有限单元与SPH粒子进行耦合。当FEM单元与SPH粒子接触时，其接触条件为龙格-库塔条件

(4)

式中：g为间隙函数；t为接触力。

SPH和FEM耦合首先是进行接触搜索，寻找临近单元，LS-DYNA中的自动接触算法是基于段的接触搜索方式，在计算过程中，SPH被定义为从节点，FEM被定义为主段。采用点面接触来实现两者的耦合，

它是通过罚函数算法实现将从节点的力作用到有限单元上的，罚函数的基本原理相当于在SPH和FEM之间加上法向接触弹簧限制质点穿透主面。在每一个计算时间步检查各节点是否穿透主面，若没有穿透不做处理，若穿透则引入接触力，其大小与主面刚度和穿透深度成正比，具体计算方法[25]不在此赘述。有限元单元与光滑粒子界面接触模型[26]如图2所示。耦合计算流程如图3所示。张志春等[27]借鉴了无网格接触算法的思想，提出了背景粒子搜索方式，将有限元节点视为背景粒子，被动地被SPH粒子搜索，任何位于有限元节点支持域内的SPH粒子都会对节点产生接触力。本文的接触搜索方式与背景粒子搜索方式不同，且不需要根据耦合界面附近有限元单元的尺寸调整耦合界面附近SPH粒子的光滑长度。

在求解耦合方程时，SPH采用条件稳定的跳蛙显式积分方法，FEM也采用条件稳定的中心差分法，两者耦合要求两者积分必须同步，这就要求两者在同一计算框架下每一步计算采用相同的计算步长，时间步长ΔtSPH-FEM取两者的较小值(式(5))。

ΔtSPH-FEM=min(ΔtSPH,ΔtFEM)

(5)

式中：SPH时间步长ΔtSPH=βh/c; FEM时间步长ΔtFEM≤Lmin/c,其中c为材料声速，β为时间步长比例系数，Lmin为最小单元尺寸。

图2 SPH与FEM接触

图3 SPH-FEM接触算法流程

2 拱棚洞冲击试验

2.1 试验模型描述

图4为钢筋混凝土拱棚洞足尺冲击试验模型,通过吊车提升质量为10 000 kg冲击器到2.5 m，5 m，10 m，20 m来分别获得250 kJ, 500 kJ, 1 000 kJ, 2 000 kJ的冲击能量。冲击器为一个组合体，上半部分为圆柱体，高度0.95 m，直径为1.25 m，下半部分为一个被切割后的球体，高度为0.3 m，直径为2 m(图5(a))。其中砂垫层被装在一个木箱中，高度0.9 m，跨度10 m，沿道路轴向5 m。拱形棚洞断面的几何尺寸和钢筋分布位置如图5所示，其中，D19表示直径为19 mm的钢筋。其他具体钢筋连接设计细节参考日本规范[28]。纵向方向沿道路长为6 m。

图4 拱棚洞冲击试验模型[10]

(a) 断面尺寸图

(b) 钢筋分布图

2.2 测试方法

冲击力时程曲线用一个量程为500g，采集频率为5 000 Hz，安装在冲击器内部的加速度计采集。位移计安装在拱形棚洞的拱圈内表面，测量其竖直方向位移，量程和采集频率分别为500 mm和1 000 Hz。

3 计算模型

3.1 计算单元

LS-DYNA在分析结构受冲击时的几何非线性、边界非线性、材料非线性计算上有着独特的优势。本文采用其进行耦合计算，计算模型如图6所示。被冲击砂土采用SPH粒子，被冲击区域尺寸按经验调试确定，取为2 m(x)×2 m(y)×0.6 m(z)，相邻粒子间距0.05 m；非冲击区域的砂土采用8节点完全积分六面体单元，单元尺寸0.1 m，在这里将砂土划分为两个区域主要考虑以下两点：一个是SPH粒子计算开销大，非冲击区域冲击变形小，可用有限元网格模拟；另一个是由于SPH的光滑核函数不具有Kronecher delta 张量的性质，且边界粒子本身存在缺陷，不适合作为边界。混凝土、岩土、冲击器采用6节点/8节点完全积分的实体单元。混凝土单元网格及其临近的岩土网格尺寸为0.1 m，为提高计算效率，远端岩土网格尺寸较大。钢筋采用桁架单元模拟[29]。划分网格后，共有22 919个SPH粒子，219 848个实体单元，13 249个桁架单元。

图6 拱棚洞计算模型(mm)

3.2 材料模型

表1所示为各材料的基本材料特性。各材料的本构模型如下所述。

砂：砂土材料是一种多孔介质材料，在强冲击作用下，会表现出可压缩泡沫材料的力学特性，为了准确描述砂土在冲击作用下的力学行为，利用了式(6)的抛物线描述其应力-应变曲线(图7(a))。可用MAT_CRUSHABLE_FOAM进行描述。

表1 各种材料特性

(6)

式中：σsand为应力，单位为MPa，εsand为体积应变。

钢筋：采用各向同性的弹塑性本构进行模拟，塑性硬化模量取为弹性模量的0.01倍(图7(c))，初始屈服强度为338 MPa，采用Von Mises屈服准则。可用MAT_PLASTIC_KINEMATIC进行描述。

岩体：岩体采用弹性体来模拟，可用MAT_ELASTIC描述。

冲击器：由于冲击器几乎不变形，为加快计算效率，冲击器采用刚体模拟，可用MAT_RIGID描述。

(a) 砂

(b) 混凝土

图7 主要材料本构模型

Fig.7 Main material constitutive models

3.3 模型边界

冲击器与砂土SPH粒子间定义为侵蚀接触；砂土有限单元与砂土SPH粒子间定义为点面接触；砂土有限单元与混凝土单元采用面面接触；砂土有限单元与岩土单元采用面面接触；混凝土单元与岩土采用共节点连接；钢筋单元嵌固在混凝土单元内由于钢筋为小应变，可用CONSTRAINED_LAGRANGE_IN_SOLID来定义钢筋与混凝土的约束。

岩土底部采用全约束边界，岩土侧壁可采用无反射边界条件。冲击器与砂土接触时开始计时总的冲击时间定义为200 ms。

4 结果分析

4.1 数值计算与试验对比

选取冲击能量为1 000 kJ的冲击锤冲击力时程曲线、拱中点位移时程曲线与数值模拟的结果对比，如图8所示，可看出，冲击力时程曲线包括快速上升、快速下降、稳定三个阶段，可视其为一个脉冲荷载，两者的时程曲线具有较好地吻合性。选取试验与数值模拟的冲击器冲击力峰值、拱中点位移峰值进行对比，如表2所示。

通过对比发现，数值模拟结果与试验结果吻合较好，冲击力峰值和拱中点位移峰值最大误差均没有超过10%。SPH-FEM方法能很好地模拟冲击器冲击拱棚洞的整个过程。随着冲击能量的增大，冲击力峰值和拱中点位移峰值也逐渐增大。数值模拟结果与试验有一定的误差，数值模拟的冲击力峰值整体稍偏大，拱中点位移峰值稍偏小，主要由以下几点原因造成：SPH方法忽略了粒子间的摩擦耗能；砂土的缓冲性功能与密实度有关，所用材料本构无法考虑到砂土的密实度影响。虽然双线性混凝土材料本构简便，但混凝土本构非常复杂，需用试验确定，造成拱中点位移时程后半段有一定出入，但总体趋势一致。

(a) 冲击力时程

(b) 拱中点位移时程

表2 试验结果与仿真结果对比

4.2 冲击过程现象分析

图9所示，为冲击锤冲击砂土的过程，冲击能量为1 000 kJ，当开始冲击锤与砂土接触时，砂土由于受到挤压和剪切作用，开始向四周流动(T=30 ms),此时砂土应力较小；随着冲击锤与砂土接触面积的增大，砂土被冲击锤挤向四周，在砂土表面形成一个微隆起区，内部区域形成了一个砂坑，在T=125 ms时，冲击锤位移达到最大值0.825 m，此时砂土应力最大值为12.74 MPa，对应的体积应变达到0.505，说明砂土体在冲击过程中存在着大变形；在T=200 ms时，冲击过程趋于静止，静止后冲击位移为0.804 m，卸载恢复的弹性变形很小，此时砂土体卸载导致其应力也变小。整个冲击过程与试验现象有较好的一致性。

T=0 ms

T=30 ms

T=125 ms

T=200 ms

图9 冲击器冲击砂土过程

Fig.9 The process of penetration by impactor

4.3 耗能分析

图10为冲击能量1 000 kJ的冲击工况的能量时程图。在冲击过程中，冲击器的动能和势能转化为砂垫层的内能、摩擦耗能、棚洞内能等，砂垫层耗散了很大一部分能量。滑移摩擦能主要由冲击过程中冲击器与砂垫层的摩擦所致。棚洞和岩石几乎没有耗散冲击耗能，说明塑性变形很小，砂垫层起到了很好的保护作用。从表3可看出，随着冲击能量的增大，土垫层耗能逐渐减小。这是由于冲击能量较小时，钢筋混凝土和岩石处于弹性阶段，不耗散冲击能量，随着冲击能量增大，钢筋混凝土开始出现塑性变形，耗散了一部分能量。但砂土垫层耗能仍然占总初始冲击动能的85%以上，它是一种很好的缓冲耗能材料。砂土垫层使得棚洞实现了分级耗能抗冲击，小冲击能量，棚洞不损伤，大冲击能量，棚洞不倒。