向量式有限元面內逆向轉角計算方法

王雙炎，李聰，楊玉龍，龐志成，湯晗青

?

向量式有限元面內逆向轉角計算方法

王雙炎，李聰，楊玉龍，龐志成，湯晗青

(浙江大學 建筑工程學院，浙江 杭州，310058)

為提高向量式有限元方法的計算效率，以三角形平面固體單元為例，提出一種用邊轉角法計算逆向運動的面內轉角的方法，并將其與常用的質心轉角方法進行對比，比較2種方法的計算效率；通過自編程序，分別驗證邊轉角法在計算單元剛體運動、小變形、大變形及復雜結構的變形時的準確性和計算效率。研究結果表明：與質心轉角法相比，采用邊轉角法計算面內逆向轉角具有較高的準確度和計算效率，該方法對于常見結構面內逆向轉角的計算是可行的、有效的。

向量式有限元；逆向運動；平面三角形單元；計算準確度；計算效率

向量式有限元[1-3](vector form intrinsic finite element, VFIFE)是基于向量力學和數值計算所提出的新型有限元數值計算方法，其核心思想在于將構件離散成若干個質點，通過描述各質點的運動從而對構件整體進行描述。該方法的計算特點在于通過時間積分實現對各點各步的逐步循環計算。由于不存在單元剛度矩陣和矩陣奇異等問題，又無需求解繁雜的非線性方程組，所以，相比傳統有限元方法，向量式有限元在計算特殊的非線性問題有著獨特的優勢。王震等[4-6]從單元形式入手，將向量式有限元理論推廣至薄板、薄殼以及四節點實體單元。HOU等[7]在此基礎上，推導了八節點正方體單元的向量式有限元應用。在向量式有限元非線性行為研究中，如有關斷裂[8-9]、碰撞[10?11]等方面的研究，目前也有了一定的進展。但由于目前還沒有商業的向量式有限元軟件，一般需要通過Matlab程序實現向量式分析計算，而Matlab自身的運行速度具有一定的局限性。而且越是復雜的結構，劃分的單元數越多，所需的計算總時間也越多，因此，提高計算效率就顯得尤為重要。YANG等[12]在其研究中闡述了向量式有限元的計算效率。傳統的向量式有限元一般采用質心轉角法計算逆向轉動角度。為了提高計算效率，本文以理論體系較為成熟的三角形平面固體單元[13-15]為例，采用邊轉角法計算逆向轉動角度。HOU等[7,16]雖使用了邊轉角法進行計算，但未對該方法在計算過程中的計算準確度和計算效率進行詳細討論。為此，本文作者對單元逆向運動過程中，邊角轉角法與質心轉角法2種面內轉角的計算方法的計算準確度和計算效率進行比較分析；通過算例驗證使用邊轉角法計算面內逆向轉動角度的準確性及其在計算效率方面所具有的優勢。

1 向量式有限元平面固體單元理論

1.1 點值描述

傳統有限元采用1個有限個單元的集合體來描述求解的結構，而向量式有限元則采用點值描述。這些質點通過相互連接形成單元，建立了質點與單元間的拓撲關系。質點是描述構件受力變形、幾何尺寸、空間位置、邊界條件的載體。構件在運動過程中，通過描述各質點的位置和速度向量來描述固體構件的整體運動。

每1個質點的運動都遵循牛頓第二運動定律。對于質點而言，m為質點的質量，為質點的空間位置向量，則質點的運動平衡方程可表示為

1.2 途徑單元

每1個質點的運動過程都是連續的，這些質點的運動軌跡可以表示為以時間為自變量的空間位置函數。在質點運動的全運動周期過程中，用一組時間節點(0，1，2，…，t)將質點的運動軌跡離散成若干微小時間段。在a～b時間段內，若質點的運動均滿足連續的時間函數，則稱2個時間節點間的時間段為途徑單元。

時間點和途徑單元可根據結構可能產生的行為進行設定。每1個離散質點連接形成的單元在設定的途徑單元內發生的變形均為小變形。當構件有大的幾何變形時，整個變形的過程被途徑單元離散為許多小變形運動過程的疊加，從而使得描述大變形變得簡單。

1.3 逆向運動

為了計算得到單元節點內力，需要從單元節點的整體位移中扣除剛體位移得到節點的純變形量，進而由虛功原理獲得單元節點內力。

在1個途徑單元內，平面固體單元的剛體位移由剛體平動和剛體轉動2個部分組成。圖1所示為1個三角形平面固體單元在a至時刻的途徑單元內的剛體位移變化。a時刻為單元基礎構型，此時單元節點分別為1a，2a，3a。單元在時刻運動到了圖示位置，此時單元節點為1，2和3。圖1中，為節點的剛體平移向量(=1，2，3)；2和3為經過逆向平移后的節點；為扣除剛體平移后的節點位移向量。

圖1 剛體逆向平移

在局部坐標系下，和a時刻節點的位置向量分別為和a，選取任意一點作為參考點(本文選取1號節點)，將單元平移，使得2個對應參考點位置重合，則參考點經過的逆向平移向量即為單元逆向剛體平移向量：

由此得到扣除剛體平移后的節點變形量為

式中：=2，3。

對單元進行逆向剛體轉動。假設逆向轉動角度為，則各節點的剛體轉動位移為

式中：=2，3。

1.4 節點內力計算

最后，經過正向運動，將局部坐標系下的單元節點內力轉換至整體坐標系下：

式中：=1，2，3。獲得單元節點內力后，即可由中央差分公式實現循環計算。

2 逆向運動方式探究

2.1 面內轉角的計算方法與計算精度

由本文第1節的分析可知：求解逆向運動過程中的逆向轉動角度至關重要，其計算準確度和計算效率將影響整體計算結果和計算效率。文獻[17?18]中，逆向運動過程的面內轉角估算值均采用了質心對應轉角取平均值的方法進行計算(以下簡稱質心轉角法)，如圖2所示。其中，aO和dO分別為單元在a和時刻的質心位置向量；aO與dO分別為a和時刻單元各節點相對于質心的方向向量。(=1，2，3)為節點相對于自身平面單元形心的轉動角度。

式中：=1，2，3。

在逆向轉動面內，取3個節點轉動角度的平均值作為逆向轉動轉角：

“面向傳動裝置的ABB AbilityTM 狀態監測解決方案也是本次ACW的亮點之一。”ABB中國機器人及運動控制事業部負責人李剛表示，“它是業內首個集成化服務，能夠將每臺設備的關鍵運行參數集中顯示，借助ABB的技術優勢，它還能使客戶提前了解維護需求，確保設備實現理想運行狀態。總而言之，借助ABB AbilityTM 狀態監測服務，用戶可以更好地掌握如何優化設備運行，減少宕機風險，延長設備壽命，降低成本并且增加收益。”

經逆向運動扣除了單元節點位移中的單元剛體平動和剛體轉動后，a基礎構架時刻和d時刻虛擬單元間的節點位移僅包含單元純變形和殘余轉動變形。這個殘余轉動變形來源于單元變形的不均勻分布。而向量式有限元在設定途徑單元時，保證每一個途徑單元內單元變形均為小變形，單元變形接近均勻變形，這就使得單元純變形對節點內力影響為一階量，而殘余變形的影響相對純變形的影響為高階量。

因此，可采用任意方法對面內轉動角度進行估算，只需要滿足以下基本條件：假設單元的純變形為0時(單元僅存在剛體運動)，估算方法得到的轉角應為正確的轉動角度。例如，當固體平面單元發生角度為的剛體轉動且單元變形量為零時，根據質心轉角估算法，1=2=3=，=成立，是正確的剛體轉角。

基于上述理論，本文建議使用另1種方法估算面內轉角，即邊轉角計算法(見圖2(b))。三角形平面固體單元經歷了逆向平動后，參考點1已重合。此時，需進一步扣除逆向剛體轉動。與逆向剛體平動一樣，任取包含作為逆向剛體平移的參考點的單元一邊作為參考邊，以參考邊的轉動角度作為剛體逆向轉動角度。圖2中選取的邊為節點1和3所形成的單元邊，則有：

式中：1為單元在a時刻節點1和節點3所形成的邊向量；2為單元在時刻經逆向剛體平移后，節點1和節點3所形成的方向向量。假設單元僅發生剛體運動而沒有純變形，且單元剛體轉動角度為，此時，根據邊角計算法，=成立，滿足面內轉角估算方法的基本條件。

圖2 面內逆向轉角計算方法對比

2.2 不同計算方法的計算效率

向量式有限元方法在求解式(1)的運動方程時，通常采用顯式時間積分法中的中央差分法。

無初始條件下的中央差分表達式為

有初始條件下的中央差分表達式為

式中：+1為第+1步的節點位置向量；為質點質量；和分別表示節點所受到的外力和節點內力；為時間步長；1=1/(1+/2)；2=1/(1?/2)。

除此之外，在向量式有限元中，還有一項重要的工作即估算臨界步長。在平面三角形固體單元內，材料軸向應力剛度一般較高，分析單元軸向運動時所需的時間步長較短，因此，一般以單元軸向運動為基準對計算過程中的臨界步長進行考量。將單元簡化成長度為的質點和勁度系數為的彈簧單自由度體系，材料的彈性模量為，橫截面積為，在外力作用下，構件長度變化量為，則有：

=(23)

由此可見，在材料參數相同的情況下，離散質點數越多，質點連接而成的單元網格越小，所能取得的臨界步長越小。向量式有限元作為一種在處理碰撞、斷裂、接觸等結構大變形問題時具有獨特優勢的新型數值計算方法，往往需要將構件離散成數量龐大的質點組來滿足對構件變形精確描述的要求。而在內力計算過程中，為避免內力誤差迅速積累而出現結果發散現象，尤其是在使用顯示時間積分求解運動方程的情況下，時間增量很小，迭代計算的循環步數一般很龐大，此時，計算效率的提升就顯得尤為重要。

圖3所示為不同計算方法的計算流程圖。相較于傳統的質心轉角計算方法，邊轉角計算法無需計算單元質心位置，向量計算量也明顯較少。對于1個時間步長內的1個單元而言，在面內逆向轉角計算步驟中，質心轉角算法需完成10次向量加減運算、8次向量乘法運算、3次向量點乘運算、6次向量求模運算、2次標量加減運算以及1次標量乘法運算；邊轉角法僅需完成2次向量加減運算與1次向量乘法運算。由計算機工作原理可知，邊轉角計算法可以有效提高計算機的運行速度。

3 算例驗證

從上述分析可以看出，取質心與節點的向量轉角平均值作為逆向轉動角度和以邊向量轉角為基準作為逆向轉動角度均是可行的。同時，在計算效率方面，由于邊轉角法的計算量更小，在多質點多單元的情況下更具優勢。本節通過自編Matlab程序進行算例模擬，對上述邊轉角理論的正確性進行驗證。

3.1 單元剛體運動

由表1可知：2種方法在處理單元剛體運動問題時，均具有較高的準確性。當單元旋轉10圈時，質心轉角法所用時間為2.7 s，邊轉角法所用時間為2.2 s。

3.2 正方形薄板受壓問題

正方形薄板示意圖及單元劃分如圖5所示。由圖5(a)可見：正方形薄板對角受到均布壓力荷載的作用，荷載沿垂直于平面的薄板厚度方向均勻分布，為2 N/m。正方形薄板對角線長=4 m，板厚0.1 m，彈性模量=1×1011Pa，泊松比=0。

由于結構具有對稱性，只需取1/4部分作為分析對象。由圖5(b)可見：構件的1/4被離散成6個質點，質點間互相連接形成共4個單元。將向量式有限元分析結果與傳統有限元分析結果進行對比，結果如表2所示。由表2可以看出：在小變形情況下，2種面內轉角計算方式的向量式有限元分析結果均與傳統有限元的理論值較符合。質心轉角法所用時間為2.2 s，邊轉角法所用時間為1.8 s。

圖3 不同計算方法流程圖

圖4 單元剛體轉動

圖5 正方形薄板示意圖及單元劃分

3.3 懸臂梁自由端受集中荷載作用

在Euler?Benoulli梁理論[19]中，懸臂梁端部受集中荷載作用發生劇烈變形是經典的非線性有限元算例，該算例常被用于驗證非線性算法理論的正確 性[20?23]。懸臂梁算例示意圖如圖6所示。懸臂梁長=10 m，梁高=1 m，梁厚為1 m。構件離散成544個質點，以三角形平面單元形式相互連接。材料彈性模量=10 Pa，泊松比=0。當外荷載以荷載因子/2(其中，為懸臂梁橫截面的慣性矩)來衡量時，此時懸臂梁變形為大變形。由于該算例是擬靜力問題，荷載加載形式取斜坡加載方式。同時，施加虛擬阻尼使結果收斂，取虛擬阻尼系數=10.0。使用質心轉角計算法和邊轉角計算法分別進行計算，驗證大變形情況下理論的正確性。

懸臂梁荷載施加點豎向位移和橫向位移對比分別如表3和表4。表3和表4中，ABAQUS的單元劃分方式均與算例的相同。由表3和表4可知：對于質心轉角法和邊轉角法2種不同的逆向轉動計算方法，其自由端在豎直和水平方向上的位移均與理論結果較吻合，2種方法得到的結果之間的差異也很小。即使是當荷載因子達到10這樣的高度非線性情況下，依然能取得理想結果，驗證了本文第3節提出的邊轉角計算方法的計算精確性。在計算效率方面，質點轉角法最終所用時長為878 s，而在相同情況下，邊轉角法所用時長為779 s。

表1 節點坐標對比

表2 節點位移對比

圖6 懸臂梁算例示意圖

表3 懸臂梁荷載施加點豎向位移對比

表4 懸臂梁荷載施加點水平位移對比

3.4 半圓弧結構受力大變形

圖7所示為半圓弧結構示意圖，圓弧底部兩端完全固定。圓弧外徑1=21 m，內徑2=20 m，材料彈性模量為10 Pa，泊松比為0.25。將結構分別離散成108，318和860個質點進行計算。半圓弧頂部位置施加1個豎直向上的集中荷載。由于該問題是1個擬靜力問題，荷載加載方式同樣采用斜坡加載方式。取虛擬阻尼系數=10.0。構件將在荷載作用下發生 大變形，驗證工程結構條件下邊轉角計算方法的正 確性。

為驗證除荷載施加點之外，圓弧拱所離散成的其他節點位移的正確性，分別選取荷載為1.0 N和3.5 N，觀察圓弧拱拱肋上沿部分其他質點產生合位移后連接所成的拱肋變形，如圖8所示。同時，選取荷載施加質點進行定量分析。圖9(a) 所示為圓弧結構離散為不同質點數情況下荷載施加點的位移與荷載的曲線(由于質心轉角法和邊轉角法的結果差距極小，圖8和圖9(a)僅展示了邊轉角法的結果)，并將結果與商業通用有限元軟件ABAQUS分析結果進行對比(其中，圓弧結構被離散成5 992個三角形單元)。

圖7 半圓弧結構示意圖

荷載/N：(a) 1.0; (b) 3.5

由圖9(a)可知：圓弧頂點將在荷載作用下先發生非線性位移變化；隨著荷載的增加，頂點兩邊的弧形結構逐漸由曲線變為直線，位移?荷載的變化又重新趨于線性相關；自編向量式有限元程序的分析結果與ABAQUS顯示的位移變化趨勢相吻合。隨著離散質點數的增加，頂點位移逐漸接近ABAQUS非線性有限元分析結果；且在相同的單元劃分情況下，向量式有限元計算結果也與ABAQUS計算結果相吻合。

由上述分析可以看出：在保證計算準確度的前提下，邊轉角法較質心轉角法有著更高的計算效率。由于向量式分析的特點，為了獲得更精確的結果，需要將結構離散成更多的質點以達到更精確描述結構的目的。此時，運用邊轉角法進行計算所能節省的絕對時間會更多。

圖9 半圓弧結構算例計算結果

圖9(b)所示為不同質點數下Matlab程序對于2種方法所需的平均計算時間。由圖9(b)可知：對于3種不同離散情況，邊轉角法相較于質心轉角法都能節省時間12%～15%；當結構離散成860個質點時，2種方法的絕對計算時差可達約360 s。

4 結論

1) 為提高向量式有限元的程序計算效率，基于向量式有限元的基本理論和假定，采用邊轉角法計算單元在逆向面內轉動過程中的逆向轉角。

2) 通過自編的Matlab程序算例對單元剛體運動、小變形、大變形、復雜結構4種情況進行了驗算，計算結果表明邊轉角法在保證較高計算準確度的前提下，相較于傳統的質心轉角法計算效率更高。

[1] TING E C, SHIH C, WANG Y K. Fundamentals of a vector form intrinsic finite element: Part 1. Basic procedure and a plane frame element[J]. Journal of Mechanics, 2004, 20(2): 113?122.

[2] TING E C, SHIH C, WANG Y K. Fundamentals of a vector form intrinsic finite element: Part II. Plane solid elements[J]. Journal of Mechanics. 2004, 20(2): 123?132.

[3] SHIH C, WANG Y K, TING E C. Fundamentals of a vector form intrinsic finite element: Part III. Convected material frame and examples[J]. Journal of Mechanics. 2004, 20(2): 133?143.

[4] 王震, 趙陽, 胡可. 基于向量式有限元的三角形薄板單元[J]. 工程力學, 2014, 31(1): 37?45. WANG Zhen, ZHAO Yang, HU Ke. Triangular plate element based on vectored finite element[J]. Engineering Mechanics, 2014, 31(1): 37?45.

[5] 王震. 向量式有限元薄殼單元的理論與應用[D]. 杭州：浙江大學建筑工程學院, 2013: 80?85. WANG Zhen. The theory and application of a vector finite element thin shell element[D]. Hangzhou: Zhejiang University. College of Civil Engineering and Architecture, 2013: 80?85.

[6] 王震, 趙陽, 楊學林. 基于向量式有限元的實體結構非線性行為分析[J]. 建筑結構學報. 2015(3): 133?140. WANG Zhen, ZHAO Yang, YANG Xuelin. Nonlinear behavior analysis of solid structure based on vector finite element method[J]. Journal of Architectura, 2015(3): 133?140.

[7] HOU Xiangying, FANG Zongde. Solid structure analysis with large deformation of eight-node hexahedral element using vector form intrinsic finite element[J]. Advances in Structural Engineering, 2018, 21(6): 852?861.

[8] DUAN Y F, WANG S M, WANG R Z, et al. Vector form intrinsic finite element based approach to simulate crack propagation[J]. Journal of Mechanics, 2017, 33(6): 797?812.

[9] HOU Xiangying, FANG Zongde, ZHANG Xijin, et al. Static contact analysis of spiral bevel gear based on modified VFIFE (vector form intrinsic finite element) method[J]. Applied Mathematical Modeling, 2018, 60: 192?207.

[10] XU Leige, LIN Mian. Numerical modeling of the configuration of a long-distance free-spanning submarine pipeline on an uneven seabed[J]. International Journal of Offshore and Polar Engineering, 2017, 27(1): 102?111.

[11] WU T Y. Dynamic nonlinear analysis of shell structures using a vector form intrinsic finite element[J]. Engineering Structures, 2013, 56: 2028?2040.

[12] YANG Yulong, CHENG J J R, ZHANG Tuqiao. Vector form intrinsic finite element method for planar multibody systems with multiple clearance joints[J]. Nonlinear Dynamics, 2016, 86(1): 421?440.

[13] TORNABENE F, FANTUZZI N, BACCIOCCHI M, et al. Mechanical behavior of damaged laminated composites plates and shells: higher-order shear deformation theories[J]. Composite Structures, 2018, 189: 304?329.

[14] 劉人懷, 薛江紅. 復合材料層合板殼非線性力學的研究進展[J]. 力學學報, 2017, 49(3): 487?506. LIU Renhuai, XUE Jianghong. Research progress in nonlinear mechanics of laminated composite plates and shells[J]. Journal of Mechanics, 2017, 49(3): 487?506

[15] 韓強.高等板殼理論[M]. 北京: 科學出版社, 2002: 216. HAN Qian. Advanced theory of plates and shells[M]. Beijing: Science Press, 2002: 216.

[16] XU Leige, LIN Mian. Analysis of buried pipelines subjected to reverse fault motion using the vector form intrinsic finite element method[J]. Soil Dynamics and Earthquake Engineering, 2017, 93: 61?83.

[17] 王震, 趙陽, 楊學林. 薄殼結構的向量式有限元屈曲行為分析[J]. 中南大學學報(自然科學版), 2016, 47(6): 2058?2064. WANG Zhen, ZHAO Yang, YANG Xuelin. Vector form intrinsic finite element method for buckling analysis of thin-shell structures[J]. Journal of Central South University of Science and Technology, 2016, 47(6): 2058?2064.

[18] 底大鈞. 三維空間柔性桿件的三角單元向量式有限元模型研究及仿真[D]. 哈爾濱：哈爾濱工程大學理學院, 2017: 24?31. DI Dajun. Research and Simulation of triangular element vector finite element model for three-dimensional flexible members[D]. Harbin：Harbin Engineering University. College of Science, 2017: 24?31.

[19] MATTIASSON K. Numerical results from large deflection beam and frame problems analyzed by means of elliptic integrals[J]. International Journal for Numerical Methods in Engineering. 1981, 17(1): 145?153.

[20] 張鵬飛. 結構破壞行為的數值模擬計算方法研究[D]. 杭州：浙江大學建筑工程學院, 2016: 84?85. ZHANG Pengfei. Study on numerical simulation method for structural failure behavior[D]. Hangzhou: Zhejiang University. College of Civil Engineering and Architecture, 2016: 84?85.

[21] 羅堯治, 楊超. 求解平面固體幾何大變形問題的有限質點法[J]. 工程力學, 2013, 30(4): 260?268. LUO Yaozhi, YANG Chao. Finite particle method for solving large deformation problems in plane solid geometry[J]. Engineering Mechanics, 2013, 30(4): 260?268.

[22] El-Abbasi N, Meguid S A. A new shell element accounting for through-thickness deformation[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 189(3): 841?862.

[23] SZE K Y, CHAN W K, PIAN T. An eight-node hybrid-stress solid-shell element for geometric non-linear analysis of elastic shells[J]. International Journal for Numerical Methods in Engineering, 2002, 55(7): 853?878.

Research on reverse rotation angle calculation method of vector form intrinsic finite element in plane

WANG Shuangyan, LI Cong, YANG Yulong, PANG Zhicheng, TANG Hanqing

(College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China)

In order to improve the calculation efficiency of vector form intrinsic finite element (VFIFE), the side-angle-method for calculation of the reverse rotation angle of the triangular plane solid element reverse motion was proposed. By comparing the side-angle-method with the commonly used centroid-angle-method, the calculation efficiencies of these two methods were studied. The accuracy and efficiency of the side-angle-method in calculating the rigid body motion, the small deformation, the large deformation and the deformation of complex structure were discussed and verified by examples using the self compiled program. The results show that compared with conventional centroid- angle-method, the side-angle-method has good computation accuracy and computational efficiency, and it is feasible and effective for calculation of the reverse rotation angle of common structures.

vector form intrinsic finite element; reverse motion; planar triangular element; calculation accuracy; calculation efficiency

TU33+9

A

1672?7207(2019)05?1135?09

10.11817/j.issn.1672?7207.2019.05.017

2018?06?19；

2018?08?19

國家水體污染控制與治理科技重大專項(2017ZX07201004)(Project(2017ZX07201004) supported by the Major Science and Technology Program for Water Pollution and Treatment)

楊玉龍，博士，講師，從事市政工程安全分析研究；E-mail: yulongy@zju.edu.cn

(編輯 伍錦花)