基于空间位形的在役索膜结构有限元模型修正与安全评估

2022-05-05丁一凡刘宇飞樊健生刘家豪

工程力学 2022年5期

丁一凡，刘宇飞,2，樊健生,2，刘家豪

(1. 清华大学土木工程系，北京 100084；2. 清华大学土木工程安全与耐久教育部重点实验室，北京 100084)

索膜结构受力高效、造型优美，集建筑学、结构力学、材料科学、计算机技术等于一体，长期以来得到设计人员与建设单位的关注与青睐[1]。索膜结构设计理论与建造技术的飞速发展应用相比[2-3]，在役索膜结构的检测、检查、评估、鉴定技术相对落后，难以为此类工程结构的安全运营提供长期有效的技术保障。虽然部分国内外专家学者在这些方面进行了针对性的研究，但迄今为止，在役索膜结构的损伤识别与性能评估尚有诸多问题需要解决。

在役索膜结构的安全性评估，无法直接参照国内现行的结构可靠性鉴定标准。根据《民用建筑可靠性鉴定标准》(GB 50292-2015)[4]、《工业建筑可靠性鉴定标准》(GB 50144-2019)[5]中的规定，钢构件的承载力、变形、缺陷、锈蚀均对安全性等级评定均有影响；《高耸与复杂钢结构检测与鉴定标准》(GB 51008-2016)[6]规定，大跨度及空间钢结构的安全性鉴定应按结构整体性和结构承载安全性两个项目分别评定等级，上述标准并无针对索膜结构的专门规定。从结构特点上讲，索膜结构属于柔性结构，施工阶段的初始张拉应力大小与索力分布决定了结构的刚度和形状，在结构服役阶段，索、膜的损伤与张拉应力的损失同样会显著影响结构位形。因此，服役阶段的索膜结构不能直接使用设计模型进行结构性能分析与安全评价，需要全面、准确地测量结构的实际位形以及索力、损伤等指标，采用基于实测数据的模型修正方法进行结构性能分析。然而，工程实践与文献资料中未见成熟可靠的在役索膜结构模型修正技术与方法。

本文提出基于空间位形参数的在役索膜结构有限元模型修正方法以实现结构的性能评估。采用现场三维扫描获取结构的真实空间位形信息，以模型中拉索的初始预应力作为优化参数、以结构关键节点坐标的匹配作为模型修正目标函数，设计ABAQUS 与MATLAB 联合仿真优化程序，通过全局搜索与局部优化实现精细有限元模型的迭代修正。修正后模型与实际结构受力状况更为接近，索力的实测与对比结果表明，有限元模型的索力由修正前最大偏差10%～30%降低至10%以下。进一步，利用修正后有限元模型对在役结构进行状态评估并得出可靠结论，为实际工程中同类结构的安全性检测鉴定、损伤识别与性能评估以及相关研究提供参考。

1 研究现状

1.1 空间位形参数获取

在役索膜结构的位形与结构的内部应力状态具有直接联系，因此变形检测是对索膜结构进行评估的重要项目，目前获取索膜结构三维信息方法与其他空间结构并无很大差别，以下几种方法均可用于获取索膜结构的三维信息：

1) 常规大地测量法。利用常规的大地测量仪器测量方向、角度、边长、高差等技术来测定变形或三维信息，该方法灵活度、精度较高，但效率较低。此外，由于构件本身具有尺寸参数，该方法选取的目标点一般是构件表面的特征点，对于大型空间结构来说，有时通过选取的特征点不能准确地表达构件的空间位置信息。

2) 三维激光扫描。三维激光扫描能大范围、高精度、高分辨率地，以非接触的方式快速获取目标表面每个采样点的三维坐标数据。目前三维激光扫描被广泛地用于工程中的变形测量[7]、三维建模[8]等领域。该方法具有较高的采样率及较高的精度，足够用于工程测量分析。

3) 数字图像法多视角几何三维重建。采集目标图像信息，可使用基于数字图像的多视角几何三维重建方法，获取结构三维点云模型与关键节点空间坐标。现场操作相对简单、设备要求低，在测绘工作与三维空间建模[9]中应用广泛。当前各种算法的实现，逐渐提高了多视角几何重建法的精度[10]。

对于体量较大的索膜结构，采用数字图像法多视角几何三维重建的后期处理工作量大、计算需求高，因此，本文采用三维激光扫描方法直接获取索膜结构的空间信息。

1.2 有限元模型修正方法

有限元模型修正已成为在结构工程各个领域中改进数值模型的一种广泛使用的方法。如今，修正、优化后的有限元模型，已成为评估整个结构使用寿命期间结构的承载能力、安全性和可维修性的必不可少的工具。对在役索膜结构进行安全性评估的一个关键点在于建立符合结构当前服役性态的精确有限元模型，以便开展准确的承载力计算。

目前，对在役索膜结构进行损伤识别或安全性分析，部分研究以及工程实践仍是基于原始设计数据建立的结构有限元模型[11]。由于索膜结构本身的特点，实际建成结构的空间位形可能与设计模型存在差异，尤其是对于服役时间较长的在役索膜结构，采用设计模型可能会出现很大误差。对于类似索膜结构的空间预应力结构有限元模型的修正方法，既有的研究较少，而且难以对修正后模型进行评估[12]。因此，需要建立适合索膜结构的有限元模型修正方法。

有限元模型修正一般基于结构的静力或动力试验数据，通过调整模型的参数(构件的刚度、材料密度等)使得模拟结果与试验结果基本吻合，同时使得有限元模型能够更精确地反映结构的实际特性。修正方法多是采用迭代的方法对模型中的相关参数进行调整，优化的一般流程如图1 所示。迭代方法是将模型修正问题转化为优化问题，其目标函数一般是根据数值模拟结果与试验结果的差异来定义的，如：

图1 有限元模型优化流程Fig. 1 The process of finite element model updating

目前，对于有限元模型的修正方法的研究很普遍，也有较多的优秀成果，从修正方式到内部算法，基本上形成了较为完善的系统，但对于具体的结构形式，需要更加细致地调整修正方法。从修正方式来看，一般基于静力学或动力学数据进行优化，追求优化效率与准确度。翁顺、朱宏平[13]提出了基于子结构有限元模型修正方法，可大幅提高大型复杂结构的有限元模型的优化效率。Wendy 等[14]基于实测的固有频率对斜拉桥有限元模型进行修正，提高了有限元模型的分析效果。Hendrik 等[15]引入了多响应目标函数，根据实测数据优化有限元模型，显著提高了有限元模型的质量。沈雁彬等[12]提出一种基于索力敏感度分析的有限元模型修正方法，利用有限元软件计算得到索力敏感度矩阵信息，将其代入结构零状态和初始态循环迭代的过程，最终得到能够准确反映结构实际受力状态的数值模型。从优化算法来看，非线性优化算法的使用更为普遍[16-17]，同时，新算法的尝试与改进也在进行。Javier 等[18]提出基于两种优化算法(和声搜索和有效集算法)及人工神经网络的有限元模型修正方法，能够显著降低模拟计算时间。

对于索膜结构，结构位形与结构内力联系紧密，因此对索膜结构有限元模型的优化，也可以认为是通过对相关参数的优化调整，使得结构位形与结构内力的关系与实测结果相吻合。基于这一点，可以通过试验得到的结构位形信息，对模型进行修正，并以结构内力作为参照对优化结果进行评估；反之亦可用实测结构内力进行模型修正。

2 有限元模型修正方法与实现

本文提出的索膜结构有限元模型修正方法主要是基于实测的结构位形指标，选择结构中拉索的预应力作为优化参数，因为造成索膜结构位形发生变化的主要是荷载、施工时拉索预应力大小以及结构运营时的预应力损失。其中，荷载造成的影响可以通过在有限元分析过程中根据试验时的实际情况进行合理的假定来减小或消除，如根据试验现场的风力情况、气温等环境状况，对有限元分析时的风荷载等条件进行相应的设置；将拉索预应力作为优化参数则可认为，施工时施加预应力及之后的预应力损失同时考虑，优化结果可视为等效的预应力值。

2.1 有限元模型修正方法

如图2 所示，索膜结构由施工至运营可以简化为3 个阶段：初始状态为施加预应力前结构的状态，仅存在约束条件，结构内力为0；施加预应力后，结构由于受牵拉位形发生变化，结构内力改变；在外荷载(重力荷载、风荷载等)作用下，结构位形与结构内力再次改变。图中，T为各索力组成的向量，X为结构特征点的位形指标组成的向量。为了清楚表明结构位形变化，以虚线为等高线作为参照。

图2 索膜结构受力阶段Fig. 2 Loading stage of cable-membrane structure

基于对三阶段过程的分析，可以建立相应的有限元模型修正方法，如下所述。

首先选取结构中各拉索的预应力T作为优化参数，在确定结构施加的外荷载之后，便可求解得到结构特征点的位形指标，因此X可看做是关于T的函数，基于结构特征点的位形指标X构建形如式(1)的目标函数：

在确定优化参数与目标函数值后，需要针对优化问题的特点选取合适的优化算法。首先，索膜结构一般规模较大，拉索数量以及可选取的位形特征点较多；另外，结构内力与结构位形之间具有复杂的关系。可以确定该问题为多元非线性优化问题，对于类似的工程优化问题的求解方法有多种，如牛顿-高斯迭代法[19]、序列二次规划(SQP)[20]等算法。本文提出的优化方法基于SQP算法，该算法基本原理是将复杂的非线性问题转化为较简单的二次规划问题，是目前公认的求解约束非线性优化问题最有效的方法之一。针对以上优化问题，其优化原理如下所述。

而本优化问题并没有显式的优化函数用以求得梯度向量与海森矩阵，因此需要用到有限差分法(式(6))，进行近似求解。

最终，当满足给定的精度要求时，输出当前迭代得到的优化参数作为最优解。

以上方法是进行优化时的迭代过程，但有时并不能直接通过迭代过程得到最优解，往往需要使用其他算法进行预处理。对于优化参数较多的有限元模型，由于可能存在较多的局部最优值，并不能直接使用上述原理进行迭代，否则可能出现不收敛或求解结果为局部最优值的情况。对于类似的工程优化问题，一般需要预先通过全局搜索[21](包括粒子群算法、遗传算法等)的方法得到若干起点或缩小优化参数范围，之后进行局部最优化。

2.2 有限元模型修正程序实现

基于上述原理与方法，本文作者实现了ABAQUS 与MATLAB 联合仿真优化程序，基于索膜结构的空间位形指标对结构有限元模型进行参数修正，程序结构如图3 所示。该程序主要通过MATLAB 脚本实现，使用MATLAB 中的优化工具箱(optimization toolbox)进行全局搜索与优化迭代过程。每次迭代的优化参数由MATLAB 写入ABAQUS 的输入文件，目标函数值通过MATLAB调用ABAQUS 程序进行数值分析并读取输出文件进而计算得到。

图3 有限元模型优化程序结构Fig. 3 The structure of finite element model updating program

3 试验验证

为了验证上述有限元模型修正方法的可行性，对北京首都国际机场南线收费大棚索膜结构进行了现场试验。该收费大棚于2008 年4 月16 日竣工并投入使用，建筑物平面呈曲壳型，索膜结构水平投影面积3275 m2，桅杆纵向间距为166.14 m，边环索锚座横向间距为53.4 m。马路两边立有4 根钢桅杆，前、后为边环索束，每根桅杆顶部交叉下拉7 根吊索，拉住边环索，膜顶中间设置1 根谷索，谷索与边环索之间由联系索联结，每根桅杆顶部有2 根下拉索与外侧的锚座节点相连保持平衡。结构所处的地理位置与现场照片如图4、图5 所示。试验当天天气状况良好，故进行有限元分析时可忽略风荷载的影响，试验内容包括使用全站仪测量结构位形、使用激光扫描仪获取结构三维点云模型、对部分拉索进行索力测量及其他结构状况检测项目。

图4 首都国际机场南线收费大棚索膜结构地理位置Fig. 4 Location of the toll station of south line of Beijing Capital International Airport

图5 索膜结构现场照片Fig. 5 Photograph of the cable-membrane structure

3.1 量测设备

试验通过全站仪(KON-XSY-211)以及三维扫描仪(PENTAX S-3180)获取结构的位形信息。基于本结构的跨度，全站仪测距误差小于5 mm，结构特征点位置信息的测量误差在毫米级别。三维扫描仪测距、精度等参数如表1 所示，现场测量照片如图6 所示。

图6 三维扫描现场照片Fig. 6 3D scanning scene photo

表1 三维扫描仪参数Table 1 3D scanner parameters

结构索力采用振动法进行测量，使用无线智能加速度计(HCF400-A0)与数据采集仪(HGL400-A1 4G)获取拉索振动信息，进而得到索力。振动法，是基于拉索索力和自振频率之间的关系，通过测量自振频率间接得到索力，可以适当考虑钢索的松弛和弯曲[22]，误差一般小于3%。

3.2 试验测量结果

通过全站仪测量的目标点主要是结构环索与谷索上的特征点，如图7 所示，由三维激光扫描获取的点云模型如图8 所示。提取点云模型中的谷索与环索的坐标信息：首先对密集的点云模型进行5%～10%的采样，可在一定程度上去除噪声，然后针对目标点，对一定球型范围内的点进行记录，求所有记录点的重心位置作为目标点坐标，可大幅减小少量噪声点的影响。其中谷索坐标信息如表2 所示，与全站仪测量结果统一坐标系，得到的对比结果如图9 所示。通过对比可以发现，北环索与谷索两种测量结果吻合良好，之后将主要采用三维扫描得到的结构三维指标对有限元模型进行优化。

图7 全站仪测量结果示意图Fig. 7 Schematic diagram of total station measurement results

图8 结构点云模型Fig. 8 Point cloud model of the structure

图9 结构位形测量结果(高度值)Fig. 9 Structural configuration measurement results

表2 谷索坐标信息Table 2 Coordinate information of valley cables

索力检测，无论是在膜结构的建设过程中还是在其日常维护检测中，都具有举足轻重的地位。索力是否处在合理的范围内，将直接影响结构的整体受力状态和线形的平顺程度。试验使用振动法对结构索力进行测量，结果见表3。

表3 索力测量结果Table 3 Cable force measurement results

3.3 有限元模型修正

3.3.1 建立有限元模型

考虑到通用有限元软件ABAQUS 强大的非线性分析功能，有限元建模使用了ABAQUS 软件，模型如图10 所示。整个模型由梁单元、桁架单元、膜单元构成：两侧桅杆采用了梁单元进行建模，每个杆件划分为1 个单元；所有的拉索采用T3D2 的桁架单元进行模拟，每根拉索简化划分为1 个单元；膜采用了壳分类中的膜单元进行建模，与拉索的桁架单元节点协调，每个三角形膜结构使用1 个三角形单元。桅杆与接地的拉索均与地面刚接。实际结构建成时间较长，无法准确获取结构的老化损伤信息，按照设计时的材料属性进行定义。

图10 有限元模型Fig. 10 Finite element model

通过温差法对拉索施加预应力，因此在之后的优化过程中，选取温度作为优化参数，本质上还是将拉索预应力作为优化参数，结合现场试验时无风，不计风荷载，仅考虑重力荷载。

3.3.2 目标函数

由于特征点数目较多，选取部分具有代表性的点构建目标函数，同时未使用到的特征点可以用来对修正后的有限元模型进行评估。分别选取两条环索与一条谷索上各3 个特征点，如图11 所示，通过其位置信息建立目标函数，该目标函数反映了有限元计算结果与实测结果的吻合程度：

图11 特征点示意图Fig. 11 Schematic diagram of feature points

式中：XM为试验测量得到的高度坐标；XA为有限元计算得到的高度坐标，i=1,2,···,9。

3.3.3 修正流程

1) 启动优化算法。输入试验测量结果并设定优化参数初始值与可行域，设置优化算法相关参数，如收敛限值、有限差分最小步长、最大迭代次数等。

2) ABAQUS 运行。通过MATLAB 脚本将优化参数写入输入文件(INP 文件)并调用ABAQUS执行输入文件，输出结果文件(ODB 文件)。

3) 求解目标函数值。MATLAB 运行python脚本读取二进制结果文件，得到相应计算值，与实测结果进行计算得到目标函数值。

4) 优化处理。利用目标函数值进行相关运算，计算修正梯度向量、海森矩阵，得到下一个优化参数继续步骤2)循环计算，直至满足收敛限值或迭代次数限值条件。

首次局部优化完成后，程序执行全局搜索算法，基于优化参数可行域，生成一系列均匀分布的优化参数，并对其进行计算筛选出部分优质的(目标函数值相对较小)参数作为迭代起点进行局部优化，最终得到全局最优解。在这一过程中，可以通过观察全局搜索过程结果，人工选取迭代起点并缩小优化参数可行域，提高收敛速度。

3.3.4 全局搜索

由于该索膜结构较为复杂、优化参数较多，可能存在多个局部最优解，因此初始值的选取对优化结果影响重大。可通过全局搜索得到多个优质的迭代起点，同时可根据全局搜索结果手动缩小优化参数的可行域，提高迭代过程的收敛速度。全局搜索过程中，目标函数的变化如图12 所示，根据变化情况可选取部分结果较好的点进行下一步迭代计算。

图12 全局搜索过程Fig. 12 Global search process

3.3.5 局部优化

基于以上全局搜索的结果，缩小优化参数的可行域，选取优质点作为迭代起点，经过300 次左右的函数计算(主要是用于求解梯度向量、修正海森矩阵)，完成27 次迭代，最终达到预先设定的迭代精度，得到最优解，过程如图13 所示。可以发现，经过全局搜索后得到的迭代起点，在迭代过程中得到进一步优化，使得目标函数值大幅下降。目标函数的下降可以说明，前述结构中的特征点的位形信息与试验测量结果更加吻合。具体情况可以通过图14 的对比看出，修正后的变形与实测变形相差很小。在施加温差应力模拟张拉后，再对结构施加外部荷载，拉索上的特征点经历变形增大、变形减小两个过程，实测以及修正后有限元的结果中特征点1 变形较小，因此在差值较小的情况下，计算偏差百分比数值较大，未在偏差对比图中绘出。其他特征点的修正后偏差均远小于修正前偏差，有限元模型优化工作已经基本完成。

图13 迭代过程Fig. 13 Iteration process

图14 特征点变形对比Fig. 14 Feature point deformation comparison

3.3.6 模型修正结果的验证

对修正后的有限元模型进行验证，判断模型与实际结构的吻合程度，Friswell 和Mottershead[23]建议将实测结果一部分用于模型修正，另一部分用于模型评估。根据试验测量得到的索力数据，对修正后的模型进行独立验证，同时也是对本文提出的有限元模型修正方法的可行性验证。试验时共测得了8 根拉索的索力值，拉索位置如图15所示，通过实测索力值计算拉索应力，与修正前后有限元模型相应构件进行对比，结果如图16 所示。与模型修正前相比，修正后的应力计算值与实测结果偏差普遍大幅减小，其中LS21、LS22、LS31、LS42、LS51 效果明显，有限元模型的索力由修正前最大偏差10%～30%降低至10%以下。

图15 拉索位置示意图Fig. 15 Schematic diagram of cable position

图16 拉索应力对比Fig. 16 Cable stress comparison

另外，未参与模型修正的特征点也可以用来对模型进行评估，选取部分特征点对修正后的模型进行验证，如图17 所示。可以发现，与修正前相比，修正后的模型变形与实测值更加接近。从应力、应变两个方面来看，修正后模型与实际结构更加接近。因此本文提出的有限元模型修正方法可显著优化有限元模型，使其反映结构真实的服役性态，从而进一步对结构进行性能评价。

图17 特征点变形验证Fig. 17 Feature point deformation verification

3.4 结构安全性评估

基于修正后的有限元模型，对结构进行承载力分析，在最不利的荷载布置下，得到结构的应力分布如图18 所示，膜面变形云图如图19 所示，由于膜结构应力较小且会遮挡拉索的应力显示，因此应力分布图对膜单元进行了消隐处理。图中显示，应力最大的拉索处于大约1/3 跨连接边环索与谷索的位置，为1035 MPa，小于抗拉强度1770 MPa，因此该结构各拉索构件的安全性满足要求。