向量式有限元的力-位移復合控制方法

曲激婷，胡 強，宋全寶

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

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

圖1 逆向運動Fig.1 Reverse motion

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

(4)

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

圖2 正向運動Fig.2 Forward motion

2 力-位移復合控制方法

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

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

步驟2劃分結構的質點和結構單元并分別進行編號；

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

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

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

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

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

步驟8更新當前分析時刻，讀入下一時刻施加在結構上的外力，將步驟5計算的質點位置作為tn時刻質點位置；

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

步驟10后處理，輸出計算結果。

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

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

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

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

3 數(shù)值分析

算例1力—位移混合控制

一個由8結點和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。設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é)的內容，采用MATLAB編程，計算結果如圖7和8所示，D點豎向位移為向下的0.476 m，B點外加荷載為向上的7.75×106N，同時可以得到AC桿的軸力為1.141 5×107N。為了驗證結果的正確性，將D點豎向荷載P=3.0×107N與計算得到的B點豎向荷載7.75×106N同時施加到結構上，采用SAP2000進行考慮大變形的靜力非線性分析，可以得到D點豎向位移為向下的0.474 m，AC桿軸力為1.141 5×107N，與采用力-位移混合控制的向量式有限元法計算結果吻合得很好，從而驗證了本文提出的力-位移混合控制方法的正確性。

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

圖9 空間24桿星型穹頂結構Fig.9 A reticulated space structure composed of 24 members

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

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

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

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

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

圖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ōu)勢。

4 結 論

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

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

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