基于空間位形的在役索膜結構有限元模型修正與安全評估

丁一凡，劉宇飛,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) 數字圖像法多視角幾何三維重建。采集目標圖像信息，可使用基于數字圖像的多視角幾何三維重建方法，獲取結構三維點云模型與關鍵節點空間坐標?，F場操作相對簡單、設備要求低，在測繪工作與三維空間建模[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)獲取結構的位形信息?；诒窘Y構的跨度，全站儀測距誤差小于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 所示?？梢园l現，經過全局搜索后得到的迭代起點，在迭代過程中得到進一步優化，使得目標函數值大幅下降。目標函數的下降可以說明，前述結構中的特征點的位形信息與試驗測量結果更加吻合。具體情況可以通過圖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，因此該結構各拉索構件的安全性滿足要求。

圖18 結構應力分布Fig. 18 The stress distribution of the structure

圖19 膜面變形云圖(豎向撓度) /mFig. 19 Membrane surface deformation cloud map (vertical displacement)

4 結論

本文提出一種基于空間位形的在役索膜結構模型修正方法，通過現場試驗驗證了該方法的可行性，采用修正后的模型實現結構的安全評估，主要結論如下：

(1) 以結構位形構建目標函數，以拉索預應力為修正參數，通過全局搜索與局部優化可實現基于空間位形的在役柔性索膜結構的有限元模型修正。

(2) 利用現場試驗進行方法驗證，經過迭代修正，有限元模型的空間位形信息與實測數據吻合良好，更準確地反映了結構真實的服役性態。利用實測索力值進行獨立驗證，結果表明方法有效可靠且精度較高。

(3) 進行索膜結構有限元模型修正，由于優化參數較多且結構非線性特征明顯，可能會求解得到局部最優值，可以通過全局搜索篩選迭代起點、縮減優化參數可行域，進而提高優化效率與準確度。