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

?

向量式有限元面內(nèi)逆向轉(zhuǎn)角計(jì)算方法

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

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

為提高向量式有限元方法的計(jì)算效率，以三角形平面固體單元為例，提出一種用邊轉(zhuǎn)角法計(jì)算逆向運(yùn)動(dòng)的面內(nèi)轉(zhuǎn)角的方法，并將其與常用的質(zhì)心轉(zhuǎn)角方法進(jìn)行對(duì)比，比較2種方法的計(jì)算效率；通過自編程序，分別驗(yàn)證邊轉(zhuǎn)角法在計(jì)算單元?jiǎng)傮w運(yùn)動(dòng)、小變形、大變形及復(fù)雜結(jié)構(gòu)的變形時(shí)的準(zhǔn)確性和計(jì)算效率。研究結(jié)果表明：與質(zhì)心轉(zhuǎn)角法相比，采用邊轉(zhuǎn)角法計(jì)算面內(nèi)逆向轉(zhuǎn)角具有較高的準(zhǔn)確度和計(jì)算效率，該方法對(duì)于常見結(jié)構(gòu)面內(nèi)逆向轉(zhuǎn)角的計(jì)算是可行的、有效的。

向量式有限元；逆向運(yùn)動(dòng)；平面三角形單元；計(jì)算準(zhǔn)確度；計(jì)算效率

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

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

1.1 點(diǎn)值描述

傳統(tǒng)有限元采用1個(gè)有限個(gè)單元的集合體來描述求解的結(jié)構(gòu)，而向量式有限元?jiǎng)t采用點(diǎn)值描述。這些質(zhì)點(diǎn)通過相互連接形成單元，建立了質(zhì)點(diǎn)與單元間的拓?fù)潢P(guān)系。質(zhì)點(diǎn)是描述構(gòu)件受力變形、幾何尺寸、空間位置、邊界條件的載體。構(gòu)件在運(yùn)動(dòng)過程中，通過描述各質(zhì)點(diǎn)的位置和速度向量來描述固體構(gòu)件的整體運(yùn)動(dòng)。

每1個(gè)質(zhì)點(diǎn)的運(yùn)動(dòng)都遵循牛頓第二運(yùn)動(dòng)定律。對(duì)于質(zhì)點(diǎn)而言，m為質(zhì)點(diǎn)的質(zhì)量，為質(zhì)點(diǎn)的空間位置向量，則質(zhì)點(diǎn)的運(yùn)動(dòng)平衡方程可表示為

1.2 途徑單元

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

時(shí)間點(diǎn)和途徑單元可根據(jù)結(jié)構(gòu)可能產(chǎn)生的行為進(jìn)行設(shè)定。每1個(gè)離散質(zhì)點(diǎn)連接形成的單元在設(shè)定的途徑單元內(nèi)發(fā)生的變形均為小變形。當(dāng)構(gòu)件有大的幾何變形時(shí)，整個(gè)變形的過程被途徑單元離散為許多小變形運(yùn)動(dòng)過程的疊加，從而使得描述大變形變得簡(jiǎn)單。

1.3 逆向運(yùn)動(dòng)

為了計(jì)算得到單元節(jié)點(diǎn)內(nèi)力，需要從單元節(jié)點(diǎn)的整體位移中扣除剛體位移得到節(jié)點(diǎn)的純變形量，進(jìn)而由虛功原理獲得單元節(jié)點(diǎn)內(nèi)力。

在1個(gè)途徑單元內(nèi)，平面固體單元的剛體位移由剛體平動(dòng)和剛體轉(zhuǎn)動(dòng)2個(gè)部分組成。圖1所示為1個(gè)三角形平面固體單元在a至?xí)r刻的途徑單元內(nèi)的剛體位移變化。a時(shí)刻為單元基礎(chǔ)構(gòu)型，此時(shí)單元節(jié)點(diǎn)分別為1a，2a，3a。單元在時(shí)刻運(yùn)動(dòng)到了圖示位置，此時(shí)單元節(jié)點(diǎn)為1，2和3。圖1中，為節(jié)點(diǎn)的剛體平移向量(=1，2，3)；2和3為經(jīng)過逆向平移后的節(jié)點(diǎn)；為扣除剛體平移后的節(jié)點(diǎn)位移向量。

圖1 剛體逆向平移

在局部坐標(biāo)系下，和a時(shí)刻節(jié)點(diǎn)的位置向量分別為和a，選取任意一點(diǎn)作為參考點(diǎn)(本文選取1號(hào)節(jié)點(diǎn))，將單元平移，使得2個(gè)對(duì)應(yīng)參考點(diǎn)位置重合，則參考點(diǎn)經(jīng)過的逆向平移向量即為單元逆向剛體平移向量：

由此得到扣除剛體平移后的節(jié)點(diǎn)變形量為

式中：=2，3。

對(duì)單元進(jìn)行逆向剛體轉(zhuǎn)動(dòng)。假設(shè)逆向轉(zhuǎn)動(dòng)角度為，則各節(jié)點(diǎn)的剛體轉(zhuǎn)動(dòng)位移為

式中：=2，3。

1.4 節(jié)點(diǎn)內(nèi)力計(jì)算

最后，經(jīng)過正向運(yùn)動(dòng)，將局部坐標(biāo)系下的單元節(jié)點(diǎn)內(nèi)力轉(zhuǎn)換至整體坐標(biāo)系下：

式中：=1，2，3。獲得單元節(jié)點(diǎn)內(nèi)力后，即可由中央差分公式實(shí)現(xiàn)循環(huán)計(jì)算。

2 逆向運(yùn)動(dòng)方式探究

2.1 面內(nèi)轉(zhuǎn)角的計(jì)算方法與計(jì)算精度

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

式中：=1，2，3。

在逆向轉(zhuǎn)動(dòng)面內(nèi)，取3個(gè)節(jié)點(diǎn)轉(zhuǎn)動(dòng)角度的平均值作為逆向轉(zhuǎn)動(dòng)轉(zhuǎn)角：

“面向傳動(dòng)裝置的ABB AbilityTM 狀態(tài)監(jiān)測(cè)解決方案也是本次ACW的亮點(diǎn)之一?！盇BB中國機(jī)器人及運(yùn)動(dòng)控制事業(yè)部負(fù)責(zé)人李剛表示，“它是業(yè)內(nèi)首個(gè)集成化服務(wù)，能夠?qū)⒚颗_(tái)設(shè)備的關(guān)鍵運(yùn)行參數(shù)集中顯示，借助ABB的技術(shù)優(yōu)勢(shì)，它還能使客戶提前了解維護(hù)需求，確保設(shè)備實(shí)現(xiàn)理想運(yùn)行狀態(tài)?？偠灾柚鶤BB AbilityTM 狀態(tài)監(jiān)測(cè)服務(wù)，用戶可以更好地掌握如何優(yōu)化設(shè)備運(yùn)行，減少宕機(jī)風(fēng)險(xiǎn)，延長(zhǎng)設(shè)備壽命，降低成本并且增加收益?！?/p>

經(jīng)逆向運(yùn)動(dòng)扣除了單元節(jié)點(diǎn)位移中的單元?jiǎng)傮w平動(dòng)和剛體轉(zhuǎn)動(dòng)后，a基礎(chǔ)構(gòu)架時(shí)刻和d時(shí)刻虛擬單元間的節(jié)點(diǎn)位移僅包含單元純變形和殘余轉(zhuǎn)動(dòng)變形。這個(gè)殘余轉(zhuǎn)動(dòng)變形來源于單元變形的不均勻分布。而向量式有限元在設(shè)定途徑單元時(shí)，保證每一個(gè)途徑單元內(nèi)單元變形均為小變形，單元變形接近均勻變形，這就使得單元純變形對(duì)節(jié)點(diǎn)內(nèi)力影響為一階量，而殘余變形的影響相對(duì)純變形的影響為高階量。

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

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

式中：1為單元在a時(shí)刻節(jié)點(diǎn)1和節(jié)點(diǎn)3所形成的邊向量；2為單元在時(shí)刻經(jīng)逆向剛體平移后，節(jié)點(diǎn)1和節(jié)點(diǎn)3所形成的方向向量。假設(shè)單元僅發(fā)生剛體運(yùn)動(dòng)而沒有純變形，且單元?jiǎng)傮w轉(zhuǎn)動(dòng)角度為，此時(shí)，根據(jù)邊角計(jì)算法，=成立，滿足面內(nèi)轉(zhuǎn)角估算方法的基本條件。

圖2 面內(nèi)逆向轉(zhuǎn)角計(jì)算方法對(duì)比

2.2 不同計(jì)算方法的計(jì)算效率

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

無初始條件下的中央差分表達(dá)式為

有初始條件下的中央差分表達(dá)式為

式中：+1為第+1步的節(jié)點(diǎn)位置向量；為質(zhì)點(diǎn)質(zhì)量；和分別表示節(jié)點(diǎn)所受到的外力和節(jié)點(diǎn)內(nèi)力；為時(shí)間步長(zhǎng)；1=1/(1+/2)；2=1/(1?/2)。

除此之外，在向量式有限元中，還有一項(xiàng)重要的工作即估算臨界步長(zhǎng)。在平面三角形固體單元內(nèi)，材料軸向應(yīng)力剛度一般較高，分析單元軸向運(yùn)動(dòng)時(shí)所需的時(shí)間步長(zhǎng)較短，因此，一般以單元軸向運(yùn)動(dòng)為基準(zhǔn)對(duì)計(jì)算過程中的臨界步長(zhǎng)進(jìn)行考量。將單元簡(jiǎn)化成長(zhǎng)度為的質(zhì)點(diǎn)和勁度系數(shù)為的彈簧單自由度體系，材料的彈性模量為，橫截面積為，在外力作用下，構(gòu)件長(zhǎng)度變化量為，則有：

=(23)

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

圖3所示為不同計(jì)算方法的計(jì)算流程圖。相較于傳統(tǒng)的質(zhì)心轉(zhuǎn)角計(jì)算方法，邊轉(zhuǎn)角計(jì)算法無需計(jì)算單元質(zhì)心位置，向量計(jì)算量也明顯較少。對(duì)于1個(gè)時(shí)間步長(zhǎng)內(nèi)的1個(gè)單元而言，在面內(nèi)逆向轉(zhuǎn)角計(jì)算步驟中，質(zhì)心轉(zhuǎn)角算法需完成10次向量加減運(yùn)算、8次向量乘法運(yùn)算、3次向量點(diǎn)乘運(yùn)算、6次向量求模運(yùn)算、2次標(biāo)量加減運(yùn)算以及1次標(biāo)量乘法運(yùn)算；邊轉(zhuǎn)角法僅需完成2次向量加減運(yùn)算與1次向量乘法運(yùn)算。由計(jì)算機(jī)工作原理可知，邊轉(zhuǎn)角計(jì)算法可以有效提高計(jì)算機(jī)的運(yùn)行速度。

3 算例驗(yàn)證

從上述分析可以看出，取質(zhì)心與節(jié)點(diǎn)的向量轉(zhuǎn)角平均值作為逆向轉(zhuǎn)動(dòng)角度和以邊向量轉(zhuǎn)角為基準(zhǔn)作為逆向轉(zhuǎn)動(dòng)角度均是可行的。同時(shí)，在計(jì)算效率方面，由于邊轉(zhuǎn)角法的計(jì)算量更小，在多質(zhì)點(diǎn)多單元的情況下更具優(yōu)勢(shì)。本節(jié)通過自編Matlab程序進(jìn)行算例模擬，對(duì)上述邊轉(zhuǎn)角理論的正確性進(jìn)行驗(yàn)證。

3.1 單元?jiǎng)傮w運(yùn)動(dòng)

由表1可知：2種方法在處理單元?jiǎng)傮w運(yùn)動(dòng)問題時(shí)，均具有較高的準(zhǔn)確性。當(dāng)單元旋轉(zhuǎn)10圈時(shí)，質(zhì)心轉(zhuǎn)角法所用時(shí)間為2.7 s，邊轉(zhuǎn)角法所用時(shí)間為2.2 s。

3.2 正方形薄板受壓?jiǎn)栴}

正方形薄板示意圖及單元?jiǎng)澐秩鐖D5所示。由圖5(a)可見：正方形薄板對(duì)角受到均布?jí)毫奢d的作用，荷載沿垂直于平面的薄板厚度方向均勻分布，為2 N/m。正方形薄板對(duì)角線長(zhǎng)=4 m，板厚0.1 m，彈性模量=1×1011Pa，泊松比=0。

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

圖3 不同計(jì)算方法流程圖

圖4 單元?jiǎng)傮w轉(zhuǎn)動(dòng)

圖5 正方形薄板示意圖及單元?jiǎng)澐?/p>

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

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

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

表1 節(jié)點(diǎn)坐標(biāo)對(duì)比

表2 節(jié)點(diǎn)位移對(duì)比

圖6 懸臂梁算例示意圖

表3 懸臂梁荷載施加點(diǎn)豎向位移對(duì)比

表4 懸臂梁荷載施加點(diǎn)水平位移對(duì)比

3.4 半圓弧結(jié)構(gòu)受力大變形

圖7所示為半圓弧結(jié)構(gòu)示意圖，圓弧底部?jī)啥送耆潭?。圓弧外徑1=21 m，內(nèi)徑2=20 m，材料彈性模量為10 Pa，泊松比為0.25。將結(jié)構(gòu)分別離散成108，318和860個(gè)質(zhì)點(diǎn)進(jìn)行計(jì)算。半圓弧頂部位置施加1個(gè)豎直向上的集中荷載。由于該問題是1個(gè)擬靜力問題，荷載加載方式同樣采用斜坡加載方式。取虛擬阻尼系數(shù)=10.0。構(gòu)件將在荷載作用下發(fā)生 大變形，驗(yàn)證工程結(jié)構(gòu)條件下邊轉(zhuǎn)角計(jì)算方法的正 確性。

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

圖7 半圓弧結(jié)構(gòu)示意圖

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

由圖9(a)可知：圓弧頂點(diǎn)將在荷載作用下先發(fā)生非線性位移變化；隨著荷載的增加，頂點(diǎn)兩邊的弧形結(jié)構(gòu)逐漸由曲線變?yōu)橹本€，位移?荷載的變化又重新趨于線性相關(guān)；自編向量式有限元程序的分析結(jié)果與ABAQUS顯示的位移變化趨勢(shì)相吻合。隨著離散質(zhì)點(diǎn)數(shù)的增加，頂點(diǎn)位移逐漸接近ABAQUS非線性有限元分析結(jié)果；且在相同的單元?jiǎng)澐智闆r下，向量式有限元計(jì)算結(jié)果也與ABAQUS計(jì)算結(jié)果相吻合。

由上述分析可以看出：在保證計(jì)算準(zhǔn)確度的前提下，邊轉(zhuǎn)角法較質(zhì)心轉(zhuǎn)角法有著更高的計(jì)算效率。由于向量式分析的特點(diǎn)，為了獲得更精確的結(jié)果，需要將結(jié)構(gòu)離散成更多的質(zhì)點(diǎn)以達(dá)到更精確描述結(jié)構(gòu)的目的。此時(shí)，運(yùn)用邊轉(zhuǎn)角法進(jìn)行計(jì)算所能節(jié)省的絕對(duì)時(shí)間會(huì)更多。

圖9 半圓弧結(jié)構(gòu)算例計(jì)算結(jié)果

圖9(b)所示為不同質(zhì)點(diǎn)數(shù)下Matlab程序?qū)τ?種方法所需的平均計(jì)算時(shí)間。由圖9(b)可知：對(duì)于3種不同離散情況，邊轉(zhuǎn)角法相較于質(zhì)心轉(zhuǎn)角法都能節(jié)省時(shí)間12%～15%；當(dāng)結(jié)構(gòu)離散成860個(gè)質(zhì)點(diǎn)時(shí)，2種方法的絕對(duì)計(jì)算時(shí)差可達(dá)約360 s。

4 結(jié)論

1) 為提高向量式有限元的程序計(jì)算效率，基于向量式有限元的基本理論和假定，采用邊轉(zhuǎn)角法計(jì)算單元在逆向面內(nèi)轉(zhuǎn)動(dòng)過程中的逆向轉(zhuǎn)角。

2) 通過自編的Matlab程序算例對(duì)單元?jiǎng)傮w運(yùn)動(dòng)、小變形、大變形、復(fù)雜結(jié)構(gòu)4種情況進(jìn)行了驗(yàn)算，計(jì)算結(jié)果表明邊轉(zhuǎn)角法在保證較高計(jì)算準(zhǔn)確度的前提下，相較于傳統(tǒng)的質(zhì)心轉(zhuǎn)角法計(jì)算效率更高。

[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]. 工程力學(xué), 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] 王震. 向量式有限元薄殼單元的理論與應(yīng)用[D]. 杭州：浙江大學(xué)建筑工程學(xué)院, 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] 王震, 趙陽, 楊學(xué)林. 基于向量式有限元的實(shí)體結(jié)構(gòu)非線性行為分析[J]. 建筑結(jié)構(gòu)學(xué)報(bào). 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] 劉人懷, 薛江紅. 復(fù)合材料層合板殼非線性力學(xué)的研究進(jìn)展[J]. 力學(xué)學(xué)報(bào), 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] 韓強(qiáng).高等板殼理論[M]. 北京: 科學(xué)出版社, 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] 王震, 趙陽, 楊學(xué)林. 薄殼結(jié)構(gòu)的向量式有限元屈曲行為分析[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 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]. 哈爾濱：哈爾濱工程大學(xué)理學(xué)院, 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] 張鵬飛. 結(jié)構(gòu)破壞行為的數(shù)值模擬計(jì)算方法研究[D]. 杭州：浙江大學(xué)建筑工程學(xué)院, 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] 羅堯治, 楊超. 求解平面固體幾何大變形問題的有限質(zhì)點(diǎn)法[J]. 工程力學(xué), 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

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

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

(編輯 伍錦花)