曲激婷，胡 強，宋全寶

(大連理工大學(xué) 建設(shè)工程學(xué)部工程力學(xué)研究所，遼寧 大連 116024)

2004年普渡大學(xué)Ting等[1-2]從向量式力學(xué)的基本概念出發(fā)，提出一種新的數(shù)值分析方法—向量式有限元(也稱有限質(zhì)點法)。該方法將結(jié)構(gòu)離散為一群質(zhì)點，用質(zhì)點的運動描述結(jié)構(gòu)的行為，同時引入“虛擬逆向運動”和“途徑單元”的概念，在不斷更新的基礎(chǔ)架構(gòu)上計算單元內(nèi)力。由于不需要集裝剛度矩陣及迭代求解非線性方程，該方法可以有效處理結(jié)構(gòu)大變形、彈塑性等非線性行為；此外，由于在計算過程中可隨意增加單元和改變邊界條件，在分析結(jié)構(gòu)不連續(xù)行為時，與傳統(tǒng)的有限元方法相比具有很大優(yōu)勢[3]，目前已在薄膜屈曲褶皺[4]、斷裂模擬[5]、橋梁倒塌[6]等多個領(lǐng)域取得了研究成果。

力控制和位移控制是實驗中兩種常用加載方法，當(dāng)結(jié)構(gòu)剛度很大時通常采用力控制，當(dāng)剛度衰減、結(jié)構(gòu)軟化時可使用位移控制[7-8]。然而，由于實驗花費高、周期長、混合控制加載方式存在困難等原因，使得基于理論分析的數(shù)值計算方法成為研究人員進行結(jié)構(gòu)分析的重要手段。目前，絕大部分基于向量式有限元的研究中，同一結(jié)構(gòu)、同一分析過程外加激勵(包括荷載和位移)的施加方法主要有兩種，即單獨使用力控制加載或單獨使用位移控制加載，這在一些實際問題中存在一定的局限性。有時，同一實驗過程中會先后使用了兩種控制加載方法；此外，實際結(jié)構(gòu)的工作狀況通常很復(fù)雜，在地震作用下建筑或橋梁結(jié)構(gòu)部分構(gòu)件會產(chǎn)生大變形、大轉(zhuǎn)動甚至斷裂，分析構(gòu)件性能變化或模擬結(jié)構(gòu)破壞過程時，通常采用力控制法在上部結(jié)構(gòu)構(gòu)件上施加自重等均布荷載，在地震發(fā)生過程中將地震加速度時程轉(zhuǎn)化為位移時程施加在結(jié)構(gòu)支座，即同一結(jié)構(gòu)的不同部分需要同時進行力控制加載和位移控制加載。向量式有限元可以實現(xiàn)實驗中的混合控制加載方式，程序編寫簡單方便，改變邊界條件靈活。目前關(guān)于力-位移混合控制的研究還較少，因此，有必要針對向量式有限元理論進行力-位移復(fù)合控制研究，為實驗加載裝置穩(wěn)定、安全、方便的實現(xiàn)力與位移混合控制提供理論依據(jù)。

本文基于向量式有限元基本原理，提出力-位移交替控制方法和力-位移混合控制的方法(統(tǒng)稱為力-位移復(fù)合控制方法)；前者考慮兩種控制方法在各自適用的階段獨立工作，根據(jù)分析過程中結(jié)構(gòu)剛度的變化，“接力”交替采用，即實現(xiàn)力控制和位移控制時間上的組合；后者是在分析過程中部分結(jié)構(gòu)采用力控制，部分結(jié)構(gòu)采用位移控制，實現(xiàn)力控制和位移控制空間上的組合。

1 向量式有限元的基本理論

向量式有限元是結(jié)構(gòu)行為分析的新方法，創(chuàng)新性地提出了點值描述、途徑單元和虛擬逆向運動的概念。該方法在空間上將結(jié)構(gòu)離散為一群質(zhì)點，用質(zhì)點的運動來描述結(jié)構(gòu)的幾何形狀和空間位置的變化，質(zhì)點之間采用無質(zhì)量的結(jié)構(gòu)單元連接和傳遞內(nèi)力，質(zhì)點的運動滿足牛頓第二定律。結(jié)構(gòu)中質(zhì)點的運動方程[9]可寫成

(1)

將式(1)轉(zhuǎn)化為差分形式進行數(shù)值求解，即

(2)

式中：d為質(zhì)點位移；Δt為時間步長；下標(biāo)n+1，n和n-1分別為第n+1，第n和第n-1時刻的變量。通過式(2)可以在外力已知的情況下求位移，即位移控制加載的方法。

當(dāng)已知結(jié)構(gòu)運動路徑時，可選用位移控制加載法計算作用在節(jié)點上的外力[10]，即

(3)

忽略慣性項對變形的影響，則結(jié)構(gòu)的實際變形可以分解為剛體運動和純變形，內(nèi)力只與純變形相關(guān)，而與剛體運動無關(guān)。因此，在求解單元的內(nèi)力時，需要從結(jié)構(gòu)的實際變形中扣除剛體運動部分,向量式有限元通過逆向運動來實現(xiàn)這一過程。以一個桁架單元為例(見圖1)，t1時刻桁架桿單元位于A1B1位置,t2時刻單元位于A2B2位置。為了求得此過程中桿單元軸力變化量，先將單元從A2B2平移以消除剛體平動的影響，再將單元以A1為參考點，逆向旋轉(zhuǎn)Δθ到A1B4位置以消除剛體轉(zhuǎn)動的影響。

圖1 逆向運動Fig.1 Reverse motion

比較A1B1和A1B4的形狀,當(dāng)Δt取值很小時，桿件在此時間段內(nèi)長度的變化量很小，可采用材料力學(xué)公式求出軸力增量Δf。再將Δf與t1時刻的桿件軸力相加，可得變形后的單元軸力為

(4)

式中：ft1，σt1，lt1，et1分別為t1時刻桿單元的軸力、應(yīng)力、長度和方向向量，ft2和lt2為t2時刻桿單元的軸力和長度。根據(jù)靜力平衡準(zhǔn)則，可求得桿件兩端其余內(nèi)力。再經(jīng)過正向的剛體平移和Δθ的剛體正向轉(zhuǎn)動,使桿件回到t2時刻位置(見圖2)，此時得到的桿件軸力為t2時刻桿件真實的軸力。

圖2 正向運動Fig.2 Forward motion

2 力-位移復(fù)合控制方法

以力控制方法為例，向量式有限元的分析流程可以概括為以下幾個步驟：

步驟1輸入各構(gòu)件的材料、截面性質(zhì)、阻尼系數(shù)、時間步長和總分析時間等參數(shù)；

步驟2劃分結(jié)構(gòu)的質(zhì)點和結(jié)構(gòu)單元并分別進行編號；

步驟3根據(jù)材料、截面性質(zhì)和構(gòu)件尺寸計算質(zhì)點質(zhì)量；

步驟4將tn時刻的質(zhì)點內(nèi)力、外力、位移、速度初始化為零向量，利用中心差分法計算tn-1時刻質(zhì)點位置和速度；

步驟5用中心差分法計算出tn+1時刻質(zhì)點位置；

步驟6令支座處質(zhì)點的tn+1時刻位置為初始位置；

步驟7遍歷所有單元，求單元在[tn,tn+1]時間內(nèi)的純變形，計算作用在質(zhì)點上的內(nèi)力；

步驟8更新當(dāng)前分析時刻，讀入下一時刻施加在結(jié)構(gòu)上的外力，將步驟5計算的質(zhì)點位置作為tn時刻質(zhì)點位置；

步驟9判斷當(dāng)前時刻是否已達到總計算時間，若否，則重復(fù)步驟5～步驟8；若是，則進入下一步；

步驟10后處理，輸出計算結(jié)果。

采用位移控制的分析流程與采用力控制方法相似，僅需在步驟6將指定自由度上的位移指定為當(dāng)前時刻的位移，在步驟8使用式(3)計算在該自由度上施加的外力，同時將“讀入下一時刻施加在結(jié)構(gòu)上的外力”修改為“讀入下一時刻施加在結(jié)構(gòu)上的位移”，其余步驟均相同。

力-位移混合控制是指在分析某種復(fù)雜工況下結(jié)構(gòu)的內(nèi)力和變形時，部分結(jié)構(gòu)采用力控制加載，部分結(jié)構(gòu)采用位移控制加載。以圖3所示的剛架為例，假定tn時刻該結(jié)構(gòu)在A點外力Pn和B點水平位移Un的共同作用下發(fā)生變形，Un是一個按照指定規(guī)律變化的位移，Un與Pn互相獨立，結(jié)構(gòu)的變形和內(nèi)力將是這兩種外加激勵共同作用的結(jié)果。

力-位移交替控制是指進行結(jié)構(gòu)分析時部分時段使用力控制，部分時段使用位移控制，但在任一個時段，結(jié)構(gòu)上只存在一種控制方法。這一方法需要解決兩個問題：一是兩種控制切換的時機；二是兩種控制切換的方法。本文對兩種控制方法的切換時機采用人為指定的方法確定，采用力控制時，當(dāng)單位力增量引起的位移大于某一量值時，將力控制切換到位移控制；采用位移控制時，當(dāng)單位位移增量對應(yīng)的外力增量大于某一量值時，將位移控制切換為力控制；起始階段可任意指定一種控制方法。本文側(cè)重討論兩種控制方法在同一分析過程中并存的可行性，因此重點討論第二個問題—兩種控制的切換方法。按照前述的分析流程，對于兩種控制的切換采用如下方法：(a)當(dāng)判定結(jié)構(gòu)需要從力控制切換到位移控制時，在步驟8記錄當(dāng)前時刻的位移x0，并從此刻起不再讀入外力增量，從下一時刻的步驟6起，在x0的基礎(chǔ)上讀入位移增量；(b)當(dāng)判定結(jié)構(gòu)需要從位移控制切換到力控制時，在步驟8記錄當(dāng)前位移增量對應(yīng)的外力P0，并從此刻起在P0的基礎(chǔ)上讀入外力增量；從下一時刻的步驟6起，不再讀入位移增量。

圖3 力-位移混合控制Fig.3 Force-displacement mixed control method

3 數(shù)值分析

算例1力—位移混合控制

一個由8結(jié)點和13個桿件組成的平面桁架[11](見圖4)，各桿件截面面積和材料特性相同，截面面積為2.5×10-3m2，彈性模量為2.1×1011N/m2，密度為7.8×103kg/m3。采用基于向量式有限元的力—位移混合控制方法，分析該平面桁架在外力和位移同時作用下的靜力大變形。時間步長取Δt=1.0×10-4s，分析時間取0.3 s，阻尼系數(shù)μ取100。設(shè)D點受豎直向下的集中荷載P=3×107N，同時在B點施加Δ=0.5 m的豎向位移(見圖5)。荷載和位移的施加過程如圖6所示。

圖4 平面桁架Fig.4 A plane truss

圖5 外力和位移施加位置Fig.5 Position of external load and displacement

根據(jù)“2”和“3”節(jié)的內(nèi)容，采用MATLAB編程，計算結(jié)果如圖7和8所示，D點豎向位移為向下的0.476 m，B點外加荷載為向上的7.75×106N，同時可以得到AC桿的軸力為1.141 5×107N。為了驗證結(jié)果的正確性，將D點豎向荷載P=3.0×107N與計算得到的B點豎向荷載7.75×106N同時施加到結(jié)構(gòu)上，采用SAP2000進行考慮大變形的靜力非線性分析，可以得到D點豎向位移為向下的0.474 m，AC桿軸力為1.141 5×107N，與采用力-位移混合控制的向量式有限元法計算結(jié)果吻合得很好，從而驗證了本文提出的力-位移混合控制方法的正確性。

圖6 荷載和位移施加過程Fig.6 The applied process of external load and displacement

圖7 D點豎向位移Fig.7 Vertical displacement of point D

圖8 B點外加荷載Fig.8 External load at point B

算例2力-位移交替控制

采用傳統(tǒng)非線性分析中經(jīng)典算例，一個由24根桿件鉸接而成的空間桁架穹頂結(jié)構(gòu)(見圖9)，桿件彈性模量為3.03×103N/m2，截面面積為3.17×10-4mm2，密度為6.309×103kg/m3，所有桿件材料特性和截面面積相同。當(dāng)作用在結(jié)構(gòu)上的荷載達到某一臨界值時，結(jié)構(gòu)將發(fā)生屈曲。設(shè)穹頂A點作用一個集中荷載P，采用向量式有限元方法對該穹頂結(jié)構(gòu)屈曲過程進行分析，時間步長取Δt=5.0×10-4s，阻尼系數(shù)μ=300。

圖9 空間24桿星型穹頂結(jié)構(gòu)Fig.9 A reticulated space structure composed of 24 members

基于向量式有限元的基本理論建立結(jié)構(gòu)模型，分別采用力控制法和位移控制法對此算例進行分析：力控制時取P=853 N,單位力增量為0.043 N；位移控制時A點豎向位移d=50 mm,單位位移增量為2.5×10-6m。

計算A點的豎向荷載與豎向位移變化情況，得到荷載-位移關(guān)系曲線即荷載-位移平衡路徑(見圖10)。將位移控制時的計算結(jié)果與文獻[12]中的計算結(jié)果進行對比，兩者基本吻合，說明向量式有限元能夠?qū)︸讽斀Y(jié)構(gòu)的屈曲過程進行追蹤，本文方法有效且分析結(jié)果正確。

對比圖10中力控制方法和位移控制方法的結(jié)果可見，兩曲線在位移為0～8.5 mm和44.5～50 mm的區(qū)段吻合得較好，而在8.5～44.5 mm區(qū)段有很大差異。差異產(chǎn)生的原因是由于采用力控制時，A點豎向位移達到8.5 mm后，繼續(xù)增加力時微小的力增量引起了巨大的位移增量，A點及與其相連桿件突然向下凹陷，即在8.5～44.5 mm區(qū)段采用位移控制時，A點不再處于平衡狀態(tài)，該段不再是星型桁架真實的平衡路徑。圖11為采用本文方法計算得到的空間桁架穹頂結(jié)構(gòu)在A點豎向位移達到44.5 mm后的形態(tài)。此結(jié)構(gòu)的剛度不大，采用位移控制時的增量步長也已取的足夠小，僅使用位移控制法就已能夠得到完整的荷載-位移關(guān)系曲線，采用力控制法追蹤荷載-位移平衡路徑并無優(yōu)勢。

圖10 A點力和位移關(guān)系曲線Fig.10 Force-displacement curves at point A

為了說明根據(jù)向量式有限元建立的分析模型使用力和位移交替控制的可行性，在A點交替使用力控制法施加外力和位移控制法施加豎向位移以追蹤空間桁架穹頂結(jié)構(gòu)的荷載-位移平衡路徑。A點位移為0～8.5 mm時使用力控制，當(dāng)A點豎向位移達到8.5 mm時，將力控制切換成位移控制，在A點豎向位移達到40 mm時，再將位移控制切換回力控制；在此過程中分析程序并未暫?；蛑貑?，直到A點豎向位移達到50 mm分析結(jié)束，分析結(jié)果如圖12 所示。兩種控制方法切換位置處計算結(jié)果發(fā)生輕微振蕩，但由于阻尼的影響及加載過程中力增量和位移增量均較小，振蕩很快消失，荷載-位移關(guān)系曲線趨于平穩(wěn)。為減弱兩種控制方法切換位置附近的振蕩現(xiàn)象，可以將切換點的時間提前，如果將第一次切換位置由8.50 mm變?yōu)?.56 mm，得到荷載-位移關(guān)系如圖13所示。由圖13可知，荷載—位移曲線變得更加光滑，并與文獻[13]的結(jié)果十分接近。說明對向量式有限元靜力分析模型采用力-位移交替控制方法是可行的，且計算結(jié)果也較為準(zhǔn)確。

圖11 空間桁架屈曲Fig.11 Buckling of the space truss

圖12 力和位移交替控制的力-位移曲線Fig.12 Force-displacement curve by alternating

圖13 切換點在5.56 mm處時的力-位移曲線Fig.13 Force-displacement curve when the force and displacement control method switching point is 5.56 mm

兩個數(shù)值算例表明：位移加載時可自動通過第一個坡，對應(yīng)的力自動增減；力加載時位移能夠自動增減；兩者結(jié)合的最終目的是希望結(jié)合兩種加載方式的優(yōu)勢。

4 結(jié) 論

本文根據(jù)向量式有限元的基本原理，提出力-位移混合控制和力-位移交替控制方法，通過數(shù)值算例的分析，驗證了本文提出的力-位移復(fù)合控制方法的正確性和可行性，此方法可以實現(xiàn)更加復(fù)雜的分析工況，拓展了向量式有限元的應(yīng)用范圍。由于向量式有限元采用物理模式描述結(jié)構(gòu)的行為，各質(zhì)點運動方程獨立求解，無需組集剛度矩陣，無需迭代求解方程，無論是力-位移交替加載，還是力-位移混合加載，都沒有給計算帶來本質(zhì)的困難，分析程序步驟清晰、過程簡單、通用性強，簡單方便的實現(xiàn)了結(jié)構(gòu)發(fā)生大轉(zhuǎn)動、大變位時力-位移復(fù)合控制方法。

(1)在運用向量式有限元的力-位移交替控制法時，不同的控制模式切換位置對于結(jié)果的影響不大；向量式有限元的運用使得力和位移混合控制變得十分便利，在處理存在大變位的結(jié)構(gòu)復(fù)雜加載過程時具有一定優(yōu)勢。

(2)對基于向量式有限元的力-位移復(fù)合控制僅展開了初步討論，在此基礎(chǔ)上，如果能夠選取合適的指標(biāo)判斷兩種控制模式切換的時機，就能夠自動地、不需要人為干預(yù)、不需要特殊處理、且在不需迭代求解非線性方程組的情況下追蹤到更加復(fù)雜的結(jié)構(gòu)的完整荷載-位移平衡路徑。