王永亮，王建辉，张 磊

(1. 中国矿业大学(北京)，力学与建筑工程学院，北京 100083；2. 中国矿业大学(北京)，煤炭资源与安全开采国家重点实验室，北京 100083)

圆弧型曲梁作为基本构件，广泛应用于土木工程、机械工程、航空航天工程等领域中[1-3]。曲梁在工程实际中常带裂纹损伤，准确评估带裂纹损伤曲梁的动力性能是结构设计的重要考虑因素；裂纹损伤深度、数目和分布均会改变曲梁基本特性，扰动梁的频率和振型[4-6]，明确裂纹损伤对动力性能的影响，可以有效确保结构的安全使用和有针对性的加固改造。同时，利用含损伤曲梁的实际自振频率和振型可以进行裂纹识别和定位[7-8]，裂纹损伤深度、数目、位置的精准识别依赖于高精度的频率和振型解答[9-10]。为了获得梁构件的高精度自由振动解答，一些解析方法和理论模型得到发展，但仍难有效应用于变曲线线型、多裂纹损伤、各类边界条件等复杂工况[1,11-12]。

数值计算成为分析复杂结构动力性能的合理选择和重要技术，有限元法被发展和应用于求解含裂纹损伤曲梁的自振频率和振型[13-15]，但解答精度依赖于网格划分质量，解答因网格划分难免引入误差[16]。特别是裂纹损伤形成各阶振型的扰动影响，对于非均匀分布网格的有效性提出较高要求。有限元网格自适应分析方法可有效地优化网格分布，在直线梁弹性屈曲[17]、曲梁振动[18-19]、板壳振动[20]、含损伤梁振动[21]和屈曲[22]、岩体变形和断裂[23]等问题求解中展示出很好的求解效力。本文将建立圆弧形曲梁裂纹的截面损伤缺陷比拟方案，进行裂纹深度、位置、数目的模拟，引入变截面Timoshenko 梁的h型有限元网格自适应分析方法[18-19]，求解含裂纹损伤圆弧曲梁自由振动问题，得到优化的网格和满足预设误差限Tol的高精度连续阶自振频率和振型。文中给出求解含多裂纹损伤圆弧曲梁自由振动数值算例，利用得出的有无损伤频率差、振型差分析了多裂纹损伤深度、数目、分布对圆弧曲梁自振频率和振型的扰动影响。

1 圆弧曲梁裂纹损伤表征方法

考虑图1 所示含裂纹损伤平面曲梁，曲梁中性轴坐标为s，坐标系为xyz，其中x、y为曲梁平面内坐标，x沿轴线切向，y沿轴线法向，z垂直于轴线所在平面。面内振动的位移为：沿x轴位移振幅u、沿y轴位移振幅v和绕z轴的转角振幅ψz。记曲梁曲率半径为R(s)，截面剪切刚度修正系数为κ，截面面积为A(s)，对z轴惯性矩为I(s)，长度为l，梁高度为h，梁厚度为b。记材料弹性模量为E，剪切模量为G，泊松比为ν，密度为 ρ。

图1 含裂纹损伤曲梁坐标系和符号Fig. 1 Coordinate systems and symbols of cracked curved beam

本文研究曲梁中的微裂纹损伤，使得梁截面产生弱化、梁的截面属性衰减。本研究采用裂纹截面损伤缺陷比拟方法[21-22]，裂纹处截面损伤定义为：

式中，Tol为自由振动解答的预设误差限。

2 圆弧曲梁自由振动

本文研究的平面曲梁面内自由振动，该特征值问题的微分控制方程为[18-19,24]：

3 网格自适应细分加密

有限元计算存在相比当前网格解答具有更高收敛阶的超收敛点[25]，利用超收敛点结合单元拼片、高阶形函数插值技术，可以提高当前有限元解的精度，得到全域的超收敛解[21,26-27]。本文对于圆弧曲梁的自由振动问题，求得当前网格下振型(位移)的有限元解后，利用有限元后处理超收敛拼片恢复方法，得到振型的超收敛解：

式中：P为给定函数向量；a为待定系数向量。随后，利用振型解答并通过Rayleigh 商计算可以获得自振频率值[28]。引入振型超收敛解，可对当前网格下振型有限元解进行能量模形式下的误差估计[21, 25]：

利用振型误差估计，网格可以进行优化处理来降低和控制振型的误差，达到预设的解答精度。本文方法对每个有限元单元e上的振型误差进行判断，如果误差控制式(10)不满足，则表明该单元上振型解答的误差过大，需要通过进行网格优化处理，本文采用单元均匀细分加密的h型网格自适应方式来增加模型自由度、降低单元上解答的误差[21]。当前单元细分生成的新单元长度与目前误差和单元阶次相关，即利用当前误差可以估计新单元的长度：

综合以上各方法可形成如下整体计算分析方案，获得含裂纹损伤圆弧曲梁的各阶频率和振型高精度解答：

1)含裂纹损伤圆弧曲梁模型。利用裂纹损伤表征方法(式(1)～式(5))，模拟多裂纹在曲梁中的深度、数目、分布，形成含裂纹损伤圆弧曲梁模型。

2)当前网格下频率和振型有限元解。在当前有限元网格下，利用曲梁自由振动问题的有限元逆幂迭代分析方法(式(6)～式(8))，求解含裂纹损伤圆弧曲梁模型，得到频率和振型的有限元解答。

3)误差估计并加密更新有限元网格。利用网格自适应细分加密方法，对当前振型解答进行误差估计，不满足预设误差限Tol时，在裂纹损伤扰动振型区域进行网格细分加密，获得更新的加密网格(式(9)～式(11))。在更新的有限元网格下，返回步骤2)、步骤3)进行循环计算和误差估计，直到获得一套充分优化的网格和满足误差限的解答。

4 数值算例

本文方法已经编制相应的Fortran 90 语言程序代码，程序开发实施基于Microsoft Visual Studio和Intel Visual Fortran 编程软件平台。本节给出求解具有代表性的多种含裂纹圆弧曲梁自由振动数值算例，对网格自适应划分以及频率、振动解答的精确性进行讨论；对裂纹损伤深度、数目、分布等因素影响自振频率和振型扰动进行了分析，检验了本文算法的可靠性和实用性。本节所有算例均采用3 次元，初始网格采用2 个单元，给定的初始误差限为Tol=10-4。

例1. 无裂纹损伤圆弧曲梁

为检验本文方法求解无损伤圆弧曲梁的精确性和有效性，本研究对图2 所示两端固定的常截面1/4 圆弧曲梁进行求解。为便于检验计算结果的数值精度，该算例采用无量纲的纯数值计算，曲梁的基本几何与物理参数如下：

图2 无裂纹损伤1/4 圆弧曲梁模型Fig. 2 Model of a quarter of uncracked circularly curved beam

表1 所示为使用本文方法求解得到的前5 阶自振频率解答，将文献[29]中结构力学求解器采用9000 个常截面直线型单元进行求解得到的高精度解答进行对比分析，同时给出了各阶求解使用的最终单元数目和频率误差，可见本文方法在自适应网格下得到的各阶解答远小于预设误差限要求。需要指出的是，本研究对振型进行误差控制，使用振型解答并通过Rayleigh 商计算得出具有更高收敛阶的频率[28]，确保频率值亦能严格满足误差限。

表1 无裂纹损伤1/4 圆弧曲梁自振频率值Table 1 Natural frequencies of a quarter of uncracked circularly curved beam

图3 所示为使用本文方法求解得到的第1 阶、第5 阶振型解答，并在横坐标轴上标记出自适应网格的最终分布情况。为方便直观显示和对比分析，图中振型结果均进行归一化处理(令最大振型值为1)。可以看出，振型在两固定端的位移均为0 值；本文方法求解各阶振型均划分出非均匀网格，且在振型变化平缓区域使用稀疏网格、在振型变化剧烈处采用了相对细密的网格，避免了全域使用一致细密网格的冗余性。随着阶次的增加，振型复杂程度增强，第5 阶振型比第1 阶使用了更多的单元。

图3 裂纹损伤1/4 圆弧曲梁振型Fig. 3 Vibration modes of a quarter of uncracked circularly curved beam

例2. 单一裂纹损伤曲梁不同裂纹深度

为检验本文方法分析含不同裂纹深度圆弧曲梁自由振动问题的有效性，本研究对图4 所示两端简支的单裂纹损伤曲梁进行求解。该曲梁裂纹深度分别取为α=0.0、0.16、0.5，裂纹位置取为β=0.6125。该圆弧曲梁夹角为120°，裂纹位置角度为73.5°，其余的基本几何与物理参数如下：

图4 单裂纹损伤曲梁不同裂纹深度(α=0.0、0.16、0.5，β=0. 6125)模型Fig. 4 Model of curved beam with single crack in different depth cases (α=0.0, 0.16, 0.5, β =0.6125)

使用本文方法分别计算了该曲线梁在三种裂纹损伤深度工况下面内自由振动的连续前5 阶特征对，计算频率值列于表2。文献[14]结合能量方法和有限元模型、文献[6]采用物理模型对上述问题进行分析，得到频率值如表2 所示。通过对比本文方法和能量方法、物理模型求解结果，可以看出随着问题复杂程度增加(如阶次增加、裂纹深度增大)，二者因基本分析方法不同导致个别阶次的本文方法与能量方法解答误差(下划线标出)略有增大，在各类裂纹深度工况的其余阶次下均展示出良好的一致性。

表2 单裂纹损伤曲梁不同裂纹深度自振频率值Table 2 Natural frequencies of curved beam with single crack in different depth cases

图5 给出了本文方法求解裂纹深度对自振频率扰动影响结果。为方便直观显示和对比分析，图中给出频率值(左侧纵轴标识刻度)和频率差(右侧纵轴标识刻度，含裂纹损伤情况频率值与无裂纹损伤情况频率值的差值)。可以看出，频率值随阶次增加，没有出现显著性差异；频率差均为负值，可知裂纹损伤的出现降低了各阶频率值；α=0.5 时的各阶频率差降低幅度均比α=0.16 大，可见裂纹损伤深度越大，则梁截面产生弱化、梁截面属性衰减程度越大，表现为频率值的显著减低。

图5 裂纹深度对自振频率扰动影响Fig. 5 Disturbance influence of crack depth on natural frequencies

图6 所示为使用本文方法求解得到单裂纹损伤深度α=0.5 下第1 阶、第5 阶振型解答。可以看出，振型在裂纹损伤附近区域出现扰动，裂纹损伤对转动位移 ψz扰动最为明显，本文方法求解自适应划分出最终非均匀的网格，在裂纹附近区域使用了相对密集的网格来适应裂纹损伤引起振型的变化，体现了本文方法自适应划分网格对各阶振型变化的适应性。

图6 单裂纹损伤曲梁裂纹深度α=0.5 下振型Fig. 6 Vibration modes of curved beam with single crack in depth case α=0.5

为分析裂纹损伤程度对振型的扰动行为，图7所示为裂纹损伤α=0.16 时各振型与相应阶次无损伤振型的差值曲线，可以看出裂纹损伤所在局部区域对振型变化有显著影响，裂纹损伤是影响各振型分量扰动变化的主要因素。在本例工况条件下，裂纹损伤对各振型分量均出现扰动，其中转动位移 ψz扰动最大。

图7 单裂纹损伤曲梁裂纹深度α=0.16 下振型扰动Fig. 7 Vibration modes disturbance of curved beam with single crack in depth case α=0.16

图8 所示为裂纹损伤α=0.5 时各振型与相应阶次无损伤振型的差值曲线，振型差值相比α=0.16时各振型更大，裂纹损伤程度越大则振型扰动愈加剧烈。通过上述结果可以看出，振型差幅值与损伤程度相关，通过定量控制裂纹损伤量，可有效控制振型扰动。本算例检验了本文方法求解含裂纹损伤圆弧曲梁解答的精确性，以及对各类裂纹损伤深度问题的适用性。

图8 单裂纹损伤曲梁裂纹深度α=0.5 下振型扰动Fig. 8 Vibration modes disturbance of curved beam with single crack in depth case α=0.5

例3. 多裂纹损伤曲梁不同裂纹数目

含裂纹损伤曲梁的裂纹数目n是影响裂纹振动特性的又一重要因素，本研究设置如表3 所示的典型多裂纹数目(n=2、3、4)和裂纹位置工况。本例采用图4 所示的两端简支圆弧曲梁，该曲梁的几何模型和基本物理参数同式(13)。

表3 多裂纹损伤数目和位置工况Table 3 Number and location of multiple cracks damage

各工况中多裂纹为图9 所示的均匀分布形式，各裂纹间夹角分别为40°(工况Ⅰ，n=2)、30°(工况Ⅱ，n=3)、24°(工况Ⅲ，n=4)。

图9 多裂纹损伤曲梁不同裂纹数目(n=2、3、4)模型Fig. 9 Model of curved beam with multiple cracks in different number cases (n=2, 3, 4)

使用本文方法分别计算了该曲线梁在三种裂纹损伤数目工况下面内自由振动的连续前50 阶特征对，遴选典型的计算频率值列于表4。可以看出，频率值随裂纹数目的增加，除个别阶次(下划线标出)频率值略有增加，整体上呈现逐渐降低的趋势。

表4 多裂纹损伤曲梁不同裂纹数目自振频率值Table 4 Natural frequencies of curved beam with multiple cracks in different number cases

图10 给出了本文方法求解裂纹数目对前5 阶自振频率扰动影响结果，可以看出频率值随阶次增加，各裂纹数目工况没有出现显著性差异；裂纹损伤数目的增加降低了各阶频率值，一般情况下，裂纹数目越多则降低程度越大；在第3 阶、第4 阶时，出现频率差为正值的情况，即裂纹损伤数目增多反而提高频率值；在第4 阶时，出现裂纹数目为4 工况的频率差比裂纹数目为3 工况的频率差更小的现象，即裂纹增多并没有显著降低频率值。综合以上结果，可知多裂纹数目与位置同时影响频率值，增加裂纹数目整体上有降低各阶频率的趋势，但因为裂纹位置的改变在某些阶次上会出现频率值增加的现象。因此，在原多裂纹损伤位置基础上，继续增加新的裂纹或增大原有裂纹深度(如本文例2)，才能出现各阶均降低的频率值。

图10 裂纹数目对自振频率扰动影响Fig. 10 Disturbance influence of crack number on natural frequencies

图11 所示为使用本文方法求解得到多裂纹损伤曲梁不同裂纹数目下首阶振型扰动解答。可以看出，振型在各多裂纹损伤附近区域出现扰动，转动位移 ψz扰动最为显著，本文方法求解自适应划分出最终非均匀的网格，在裂纹附近区域使用了相对密集的网格来适应裂纹损伤引起振型的变化，体现了本文方法自适应划分网格对多裂纹损伤曲梁各阶振型变化的适应性。

图11 多裂纹损伤曲梁不同裂纹数目振型扰动Fig. 11 Vibration modes disturbance of curved beam with multiple cracks in different number cases

例4. 多裂纹损伤曲梁不同裂纹分布

为进一步分析多裂纹损伤分布对曲梁自由振动的影响，本例仍采用图4 所示的两端简支圆弧曲梁，该曲梁的几何模型和基本物理参数同式(13)。该曲梁设置5 条裂纹损伤，考虑裂纹损伤沿曲梁均匀分布(各裂纹位置β 为0.1667、0.3333、0.5000、0.6667、0.8333)、裂纹损伤集中于曲梁左侧集中分布(各裂纹位置β 为0.0833、0.1667、0.2500、0.3333、0.4167)2 种工况，各工况的多裂纹如图12 所示，各裂纹间夹角分别为20°(裂纹均匀分布)、10°(裂纹左侧集中分布)。

图12 多裂纹损伤曲梁不同裂纹分布模型Fig. 12 Model of curved beam with multiple cracks in different distribution cases

使用本文方法分别计算了该曲线梁在不同裂纹分布工况下面内自由振动的连续前50 阶特征对，遴选典型的计算频率值列于表5，裂纹分布对自振频率扰动影响如图13 所示。可以看出，裂纹左侧集中分布相比均匀分布，频率值在低阶(如第1 阶～第4 阶)时具有更高的数值，而在高阶时具有更低的数值，相同数目多裂纹的不同分布形式成为影响振动特性的重要因素。因此，需要同时精准检测出裂纹损伤数目和各裂纹位置，才能准确估计含裂纹损伤时的频率值。

图13 裂纹分布对自振频率扰动影响Fig. 13 Disturbance influence of crack distribution on natural frequencies

表5 多裂纹损伤曲梁不同裂纹分布自振频率值Table 5 Natural frequencies of curved beam with multiple cracks in different distribution cases

图14 所示为使用本文方法求解得到多裂纹损伤曲梁不同裂纹数目下首阶振型扰动解答。可以看出，振型在各均匀分布和集中分布的裂纹损伤附近区域出现扰动，本文方法在裂纹附近区域使用了相对密集的网格来适应裂纹损伤引起振型的变化，体现了本文方法自适应划分网格均能很好适应不同裂纹密集程度诱发的振型扰动。需要指出的是，在图14(b)所示的裂纹左侧集中分布工况下，振型在右侧区域同时出现较大幅值的振型扰动，体现了多裂纹损伤对整体振型的强扰动行为。

图14 多裂纹损伤曲梁不同裂纹分布振型扰动Fig. 14 Vibration modes disturbance of curved beam with multiple cracks in different distribution cases

5 结论

本文建立圆弧形曲梁裂纹的截面损伤缺陷比拟方案和h型有限元网格自适应分析方法，求解含裂纹损伤圆弧曲梁自由振动问题，得到优化的网格和满足预设误差限的高精度自振频率和振型解答，定量研究多裂纹损伤深度、数目、分布形式等对圆弧曲梁自振频率和振型的扰动行为。本文的主要结论如下：

(1)裂纹区域自适应网格。自适应网格算法对无损伤、含损伤曲梁分析具有良好适用性，振型在裂纹损伤附近区域出现扰动，本文自适应优化出非均匀网格，在裂纹附近区域使用了相对密集的网格来适应裂纹损伤引起振型的变化。

(2)裂纹损伤深度。裂纹损伤的出现降低了各阶频率值，损伤深度越大则降低程度越大；裂纹损伤对转动位移 ψz扰动最大，损伤程度越大越加剧扰动幅值。

(3)多裂纹损伤数目。裂纹数目与位置同时影响频率值，增加裂纹数目整体上有提高频率的趋势，但因为裂纹位置的改变在某些阶次频率上亦会降低频率值，各裂纹损伤附近区域的振型均出现扰动。

(4)多裂纹损伤分布。裂纹一侧集中分布相比均匀分布，频率值在低阶时具有更高的数值，而在高阶时具有更低的数值；振型在各均匀分布和集中分布的裂纹损伤附近区域均出现扰动，相同数目多裂纹的不同分布形式成为影响振动特性的重要因素。