考虑剪切和扭转变形的有限质点法纤维梁单元研究

林贤宏，罗尧治，唐敬哲，汪 伟，郑延丰，杨 超

(浙江大学空间结构研究中心，浙江，杭州 310058)

纤维梁单元模型为结构的弹塑性分析提供了一种高效通用且精度较高的模拟手段，在有限元中得到广泛的研究与应用。在纤维梁模型中，梁截面被划分成若干离散的纤维，且假定截面变形前后保持平面状态，通过独立分析截面上各纤维的弹塑性状态，进而实现对梁截面塑性发展过程的精细化模拟。传统纤维梁模型[1-2]中假定纤维为单轴受力状态，仅考虑纤维材料的单轴本构关系，能够较好地模拟梁单元在轴力和弯矩为主作用下的非线性行为，但在模拟以剪切变形为主的非线性行为时效果不好。Ranzo 等[3]和Marini 等[4]通过引入截面剪力与剪应变的非线性关系来考虑剪切变形的影响，模型简单且易于实现，但无法完全考虑轴力、弯矩和剪力之间的耦合作用。为了更真实地模拟各种外力作用下构件的弹塑性响应，一些学者通过引入纤维材料的多轴本构关系，推导了能考虑剪切变形[5-10]和扭转变形[11]影响的纤维梁模型。胡郑州等[12]推导了考虑剪切效应的三维纤维梁大位移增量非线性有限元UL 列式，将纤维梁单元应用到考虑大位移、大转动的几何非线性分析中，但尚未考虑材料非线性的影响。目前，纤维梁单元模型在有限元理论中得到了较好的发展与应用[13]，但是基于有限元理论发展的纤维梁单元在求解考虑大位移、大转动的几何非线性问题时需要采用TL 或UL 列式进行迭代求解，在求解结构材料非线性问题时需要进行刚度矩阵的修正与迭代，分析较为繁琐，求解代价较高。

与有限元方法不同，基于向量式力学发展的有限质点法(finite particle method，FPM)采用离散的质点来描述连续体，通过求解质点的运动方程来获得结构的响应特征，适合于求解复杂的结构非线性问题[14]。在有限质点法中，各质点运动方程的求解过程相互独立，不涉及结构刚度矩阵的集成和求逆，无需迭代求解，避免了结构刚度矩阵奇异造成的数值求解困难。该方法在机构运动[15-16]、结构动力非线性分析[17-18]、结构倒塌破坏[19]、薄壳屈曲[20]、膜结构形态分析[21-22]及结构多尺度精细化模拟[23]等问题的分析中均取得了较好的效果。目前，Tang 等[24]针对有限质点法提出了基于GPU 并行加速的计算策略，并应用于有限质点法通用计算平台的研发中，进一步提高了有限质点法的求解效率。

已有学者推导了有限质点法的梁单元，并围绕其弹塑性行为展开了研究。有限质点法在求解结构的弹塑性行为时无需进行结构刚度矩阵的修正与迭代，只需要在每个途径单元内按小变形假定根据本构关系求解单元内力，具有求解简单与计算效率高的优点。喻莹[19]推导了梁单元的质点内力求解方法，给出了利用塑性铰法分析结构的弹塑性行为的计算流程，能简单高效地实现结构宏观破坏模拟分析。但塑性铰法无法考虑塑性在截面上及沿杆长方向上的塑性分布与发展，计算精度不高。Duan 等[25]将传统纤维梁模型引入向量式有限元中，模拟了斜拉桥的倒塌破坏过程。Liu等[26]将考虑扭转变形和材料二维本构关系的纤维梁单元引入向量式有限元中，模拟了埋地管道在走滑断层附近的大变形行为。目前有限质点法中纤维梁的研究和应用还较少，现有的纤维梁单元尚无法考虑剪切变形，具有一定的局限性。

本文推导了可考虑剪切和扭转变形的有限质点法纤维梁单元。从有限质点法的基本原理出发，将考虑材料三维本构关系的纤维梁模型融入有限质点法的计算求解框架中，基于Timoshenko梁理论推导了能够考虑剪切变形和扭转变形的纤维梁单元(TM 纤维梁) 的质点内力求解公式，为处理考虑大位移、大转动的弹塑性问题提供了一种有效的解决方案，相较于塑性铰法具有更高的精度，相较于传统有限元法具有更高的计算效率。

1 有限质点法基本方法

1.1 质点运动方程

1.2 质点内力

有限质点法通过利用虚拟逆向运动排除结构刚体位移，得到单元的纯变形，将单元内力的求解限制在小变形小转动的范围内，从而可以方便地求得质点内力。以空间梁单元AB为例，建立以A为原点的单元局部变形坐标系，x轴的方向为单元的长度方向，y轴和z轴在垂直于x轴的平面内，如图1 所示。考虑该单元在时间段ta～tb内从AB运动到A′B′，质点A′和B′经过虚拟逆向运动，分别运动到A′′和B′′，从而可以求得梁单元AB的纯变形分量ud，具体的推导过程可参见文献[19]。单元纯变形分量ud可表达为：

图1 梁单元的纯变形Fig. 1 Pure deformation of beam element

由单元的纯变形可求得单元在局部变形坐标系下的相应的单元质点内力，经过坐标变换即可得到全局坐标系下的质点内力。有限质点法具有统一的计算求解框架，不同单元类型的质点内力求解方法的实现只需在相应的求解模块中进行修改，易于在程序中实现。下面将介绍TM 纤维梁单元的质点内力求解方法。

2 基于TM 纤维梁单元模型的质点内力求解

2.1 纤维梁单元模型基本原理

在纤维梁单元模型中，梁截面被离散为一定数量的纤维，几种典型截面的纤维离散方式如图2所示。纤维梁单元模型假定截面在单元变形过程中始终保持平面，假设每个纤维截面上的应变均匀分布，可取每根纤维中心点的应变作为其代表值。在得到各个纤维的应变后，结合纤维材料的本构关系可以计算得到截面各个位置的应力状态，对截面内力进行积分计算得到单元内力。纤维梁单元模型通过独立求解截面上各纤维的应力状态实现了对梁单元截面的塑性发展过程的模拟。

图2 典型截面的纤维离散Fig. 2 Fiber discretization of typical sections

有限质点法中单元的变形与内力的求解在单元局部变形坐标系中完成，下文将以梁长为L的矩形截面纤维梁单元AB(如图3 所示)为例，介绍基于Timoshenko 梁理论的有限质点法纤维梁单元模型的描述方式，并给出质点内力求解方法。

图3 纤维梁单元示意图Fig. 3 Schematic diagram of fiber beam element

2.2 单元位移描述

空间梁单元AB在局部变形坐标系下的位移场可描述为：

2.3 纤维状态的确定

在通常的空间梁单元理论中，可以只考虑截面上任一点的3 个相互独立的应力和应变分量[27]：

材料达到屈服后应变增量可分为弹性部分和塑性部分，即：

其他材料的弹塑性本构关系可根据具体情况作相应的推导与修改。

2.4 单元质点内力的求解

单元纯变形分量和与之对应的单元质点的基本内力之间的关系由虚功原理确定，它们满足以下方程：

3 算例分析

本文算法已在有限质点法计算平台[24]中编制了相应的程序模块，下面将通过具体算例来验证其在解决结构弹塑性问题中的正确性和有效性，并以单层网壳结构的动力响应分析为例比较了有限质点法平台与通用有限元软件ANSYS 的计算效率。

3.1 空间直角刚架

图4 所示的空间直角刚架结构是空间杆系结构弹塑性求解的典型算例。结构的几何和材料特性见图中参数，材料为理想弹塑性。竖向荷载F作用于BD杆的中点C，结构承受弯矩、扭矩和剪力的共同作用。分别取杆件截面边长a=0.25 m和a=4 m 这2 种情况进行计算分析。采用有限质点法求解该结构的弹塑性变形，将结构离散为41个质点(40 个单元)，将杆件截面离散为10×10的纤维截面。计算时采用力控制的加载方式，时间步长取为2×10-5s。

图4 空间直角刚架Fig. 4 Spatial right-angle frame

当a=0.25 m 时，节点B的竖向位移与荷载F之间的无量纲曲线如图5 所示，并与Park 等[29]和喻莹[19]的计算结果进行了对比。从图中可以看出，由于喻莹在分析时采用基于塑性铰法的屈服准则，结果呈现出多段线的特征；而本文所采用的有限质点法纤维梁单元能有效地模拟杆件在弯扭作用下的截面塑性的发展过程，与Park 采用有限元塑性区法的分析结果吻合较好，屈服曲线均较为光滑。

图5 空间直角刚架的荷载-位移曲线Fig. 5 Load-deflection curve of spatial right-angle frame

当a=4 m 时，节点B的竖向位移与荷载F之间的无量纲曲线如图6 所示，图中还给出了通用有限元软件ANSYS 的Beam188 单元的计算结果，并与a=0.25 m 时的结果进行了对比。结果表明：① 由于截面尺寸的增大，剪切变形的影响更为显著，该结构的无量纲荷载系数有所下降；② 在弹性阶段，有限质点法TM 纤维梁单元与Beam188 单元的模拟结果基本一致，可较好地反映剪切变形的影响；③ 在塑性阶段，Beam188 单元的模拟结果要略大于有限质点法TM 纤维梁单元的模拟结果，这是由于ANSYS 中的Beam188单元中假设剪力与剪应变为弹性关系，而有限质点法TM 纤维梁单元考虑了正应力与剪应力的耦合作用，并在计算中采用了材料的三维弹塑性本构关系。

图6 空间直角刚架的荷载-位移曲线对比Fig. 6 Load-deflection curve comparison of spatial right-angle frame

3.2 空间框架

图7 所示为柱子底端刚接的空间框架结构，其几何和材料特性见图中参数，材料取等向强化模型。水平集中荷载H作用于E点，竖向集中荷载V=1.3H作用于杆件DE的中点G。Harrison[30]对该结构进行了试验研究，Teh 等[31]采用有限元塑性区法对该结构进行了弹塑性计算分析。采用有限质点法求解该结构的弹塑性变形，将每根杆件离散为11 个质点(10 个单元)，将杆件截面离散为32 根纤维，计算时间步长取为5×10-6s。

图7 空间框架Fig. 7 Space frame

节点A的水平位移与荷载H之间的关系曲线如图8 所示，并与Harrison[30]和Teh 等[31]的结果进行了对比。从图中可以看出，本文方法能较好地模拟该结构的弹塑性变形行为，与Teh 采用有限元塑性区法分析的结果基本一致，且接近于Harrison的试验结果。

图8 空间框架的荷载-位移曲线Fig. 8 Load-deflection curve of space frame

3.3 六角空间刚架

图9 所示的六角空间刚架结构是空间杆系结构弹塑性失稳分析的典型算例，包括了几何大变形和材料非线性因素的影响。该结构的几何和材料特性见图中参数，材料为理想弹塑性。竖向集中荷载P作用于结构顶点。

图9 六角空间刚架 /mFig. 9 Framed dome

采用有限质点法模拟该结构的弹塑性失稳过程，每根杆件被均匀离散为9 个质点(8 个单元)，将杆件截面离散为10×10的纤维截面。计算时采用位移控制的加载方式跟踪结构失稳后的大变形行为，时间步长取为6×10-5s，单位位移增量Δu=6×10-6m。同时，本文还采用了通用有限元软件ANSYS 的Beam188 单元进行了模拟计算。

图10 给出了本文计算得到结构顶点竖向位移与荷载大小之间的关系曲线，并与Park 等[29]、喻莹[19]以及通用有限元软件ANSYS 的模拟结果进行了对比。由图可见，本文提出的有限质点法TM 纤维梁单元能够有效地模拟该结构的弹塑性失稳过程，与Park 采用有限元塑性区法分析的结果吻合较好，与ANSYS 软件的计算结果基本吻合，相较于喻莹采用塑性铰法分析的结果具有更高的精度。

图10 六角空间刚架顶点的荷载-位移曲线Fig. 10 Load-deflection curve of framed dome

3.4 单层球面网壳结构

图11 所示为周边铰接的K6 型单层球面网壳结构，其几何和材料特性见图中参数，材料为理想弹塑性。节点A作用有如图12 所示的冲击荷载P，分别取荷载幅值P0=10 kN和P0=100 kN进行动力响应分析。本文在有限质点法计算平台FPM 中采用TM 纤维梁单元对杆件进行模拟，将每根杆件离散为6 个单元，共有5580 个单元。将杆件截面均匀离散为16 根纤维，计算时间步长取2×10-5s，不考虑阻尼的作用。作为对比，选取ANSYS 软件对该结构进行动力响应分析，采用Beam188 单元对杆件进行模拟，将每根杆件划分为6 个单元，迭代步的最大时间步长取2×10-3s。两种计算程序均在同一台计算机上运行，其基本硬件参数如表1 所示。

表1 计算机硬件参数Table 1 Computer hardware parameters

图11 单层球面网壳结构Fig. 11 Single-layer spherical reticulated shell structure

图13 和图14 分别给出了荷载幅值P0=10 kN和荷载幅值P0=100 kN时，节点A 的竖向位移时程曲线。由图可知，本文算法与ANSYS 的模拟结果吻合较好。ANSYS 计算时采用的是基于传统有限元的隐式动力分析，处理非线性问题时需要进行迭代求解。而FPM 采用中心差分法求解运动方程，无需进行迭代，具有较高的计算效率。两种程序计算所耗费的时间如表2 所示。当P0=10 kN时，FPM 的计算耗时仅为ANSYS 的19%。当P0=10 kN时，结构的非线性特征更为显著，ANSYS需要耗费更多的时间用于迭代计算以保证结果的收敛性，FPM 的计算耗时仅为ANSYS 的10%，具有较大的优势。

图13 节点A 的竖向位移时程曲线(P0=10 kN)Fig. 13 Vertical displacement time history curve of node A (P0=10 kN)

图14 节点A 的竖向位移时程曲线(P0=100 kN)Fig. 14 Vertical displacement time history curve of node A (P0=100 kN)

表2 FPM 与ANSYS 的计算效率对比Table 2 Calculation efficiency comparison between FPM and ANSYS

4 结论

本文将纤维梁模型基本原理与有限质点法相结合，并考虑纤维材料的三维本构关系，提出了考虑剪切和扭转变形的有限质点法纤维梁单元模型，基于Timoshenko 梁理论推导了相应的单元质点内力求解公式，为处理考虑大位移、大转动的弹塑性问题提供了一种有效的解决方案。将纤维梁模型与有限质点法相结合，可以充分发挥两者的优势，在处理几何和材料双非线性问题中具有求解简单、易于实现的特点。数值算例分析表明，本文提出的纤维梁单元能较好地模拟结构弹塑性变形及结构大变形失稳等问题，且能够考虑杆件截面的塑性发展过程，相较于FPM 塑性铰法具有更高的精度，相较于传统有限元法具有更高的计算效率。