賈 延， 劉長武

（1.北方民族大學(xué) 信 息與計(jì)算科學(xué)學(xué)院，寧夏 銀 川 750021；2.四川大學(xué) 水 利水電學(xué)院，四川 成 都 610065）

動力機(jī)械在航空、冶金、核能及石油化工等許多領(lǐng)域有著廣泛的應(yīng)用，這些設(shè)備中一些關(guān)鍵構(gòu)件（如壓力容器、齒輪、汽輪機(jī)轉(zhuǎn)子、鍋爐汽包及汽輪機(jī)缸體等）工作的可靠性直接關(guān)系到整套設(shè)備的經(jīng)濟(jì)性與安全性。在這些構(gòu)件中，有一部分表現(xiàn)出軸對稱特征，如幾何形狀、物理性質(zhì)、運(yùn)動與動力特征、約束與邊界條件等均對稱于某一固定軸，如果在外荷載作用下，其應(yīng)力、應(yīng)變及位移也對稱于該軸，可按軸對稱問題來處理，為工程設(shè)計(jì)帶來了方便。因此，研究軸對稱構(gòu)件的強(qiáng)度和剛度問題，對于改進(jìn)設(shè)計(jì)、降低成本、提高其使用壽命有著重要的意義。

目前，解決工程上的軸對稱問題主要采用經(jīng)驗(yàn)設(shè)計(jì)或彈性力學(xué)中的解析法。一般采用解析方法很難得到其準(zhǔn)確的應(yīng)力分布，且解析法中常需要重復(fù)構(gòu)造應(yīng)力函數(shù)，帶有一定的偶然性和相當(dāng)大的局限性。隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值方法得到了人們的重視。在各種數(shù)值方法中，有限元方法是發(fā)展最為成熟以及應(yīng)用最為廣泛的一種，也為工程實(shí)踐提供了大量寶貴的解決方案。但是隨著現(xiàn)代工業(yè)的發(fā)展，有限元法遇到了一些挑戰(zhàn)，如動態(tài)裂紋的擴(kuò)展、材料的破壞與失效等不連續(xù)問題以及一些大變形問題，采用有限元法解決是有困難的。

無網(wǎng)格法［1－2］的出現(xiàn)突破了傳統(tǒng)有限元法基于網(wǎng)格劃分的限制，在分析動態(tài)裂紋擴(kuò)展［3］、金屬加工成型以及高速碰撞等涉及網(wǎng)格畸變、特大變形和自適應(yīng)分析［4］等方面具有突出優(yōu)勢，現(xiàn)已成為國內(nèi)外力學(xué)界的研究熱點(diǎn)。諸多學(xué)者通過多種不同途徑構(gòu)造不同的近似場函數(shù)并進(jìn)行積分，從而衍生出了多種具有各自特點(diǎn)的無網(wǎng)格法，而耦合多項(xiàng)式基的徑向點(diǎn)插值法［5－6］是其中之一。它在點(diǎn)插值方法［7］基礎(chǔ)上，引入徑向基函數(shù)（RBF）構(gòu)造近似函數(shù)，耦合的近似場函數(shù)兼并了多項(xiàng)式基函數(shù)和徑向基函數(shù)的優(yōu)良特性，確保了RBF的力矩矩陣可逆，同時(shí)其形函數(shù)滿足Kronecker－Delta函數(shù)性質(zhì)，形函數(shù)及其導(dǎo)數(shù)計(jì)算相對簡單，使得位移邊界條件的施加變得容易，從而克服了以往無網(wǎng)格法難以實(shí)現(xiàn)位移邊界條件的難點(diǎn)。

本文將耦合多項(xiàng)式基的徑向點(diǎn)插值法應(yīng)用于固體力學(xué)軸對稱問題中，得到其無網(wǎng)格離散方程，數(shù)值算例驗(yàn)證了該方法的有效性和適用性，并通過與ANSYS有限元軟件計(jì)算結(jié)果相比較，表明該方法是求解工程軸對稱問題的一種較為有效的方法。

1 耦合多項(xiàng)式基的徑向點(diǎn)插值法

1.1 形函數(shù)的構(gòu)造

在耦合多項(xiàng)式基的徑向點(diǎn)插值法中，計(jì)算域用一系列點(diǎn)來離散，每個(gè)點(diǎn)都有一定的影響域（influencing domail），在計(jì)算點(diǎn)x處，u（x）可用徑向基和多項(xiàng)式基的線性組合近似表示為：

其中，Ri（x）為徑向基函數(shù)（Radial Basis Function，簡稱RBF）；n為 R BFs的個(gè)數(shù)；ai為Ri（x）的系數(shù)；pj（x）為空間坐標(biāo)xT＝［x，y］中的單項(xiàng)式；m為多項(xiàng)式基函數(shù)的個(gè)數(shù)，若m＝0時(shí)，u（x）為單純的RBFs，否則為添加了m個(gè)多項(xiàng)式基函數(shù)的RBF。在二維問題中，一般采用線性基pT（x）＝［1，x，y］。

為避免采用多項(xiàng)式基PIM所引起的力矩矩陣奇異性問題，引入了徑向基。本文采用復(fù)合2次（MQ）函數(shù)，其形式為：

其中，ri為高斯點(diǎn)影響域中相鄰節(jié)點(diǎn)之間的距離；αc、q為形狀參數(shù)［8］；dc為與計(jì)算點(diǎn)x的局部支持域中節(jié)點(diǎn)間距有關(guān)的特征長度，該長度通常等于局部支持域中的節(jié)點(diǎn)平均間距。

在徑向基函數(shù)Ri（x）中，僅有表示計(jì)算點(diǎn)x與節(jié)點(diǎn)xi之間距離的變量為：

距離是無方向的，所以有Rk（ri）＝Ri（rk），即Rk（xi，yi）＝Ri（xk，yk）成立。

（1）式中的待定系數(shù)ai和bj可通過n個(gè)影響域內(nèi)節(jié)點(diǎn)的插值過程得到，如第k個(gè)點(diǎn)的插值方程為：

（1）式的矩陣形式為：

其中，函數(shù)值向量為

RBFs的系數(shù)向量為

多項(xiàng)式系數(shù)向量為

RBFs的力矩矩陣為

多項(xiàng)式力矩矩陣為

然而在（3）式中有n＋m個(gè)未知量，卻只有n個(gè)方程，可使用m個(gè)約束條件添加m個(gè)方程，即

聯(lián)立方程（3）和（6），則得m＋n個(gè)方程，可解m＋n個(gè)未知量，寫成矩陣形式為：

對于任意的節(jié)點(diǎn)分布，R0的逆矩陣通常是存在的，這一點(diǎn)數(shù)學(xué)家Powell和Schaback分別在1992年和1994年已證明。因?yàn)樵冢?）式中RBFs的力矩矩陣R0是對稱矩陣，故矩陣G也將是n＋m階對稱矩陣。由于建立的近似函數(shù)是局部的，所以由其建立的無網(wǎng)格方法的系數(shù)矩陣為稀疏矩陣。求解（7）式可得：

將（1）式改寫為：

將（8）式代入（9）式，最后求解得：

其中，Φ（x）＝［φ1（x），φ2（x），…，φn（x）］稱為耦合多項(xiàng)式基的徑向點(diǎn)插值法的形函數(shù)；φk（x）＝，稱作第k個(gè)節(jié)點(diǎn)的形函數(shù)i，k為矩陣G－1的元素。形函數(shù)的導(dǎo)數(shù)形式為：

1.2 離散方程的推導(dǎo)

為方便計(jì)算，軸對稱問題常采用柱坐標(biāo)（r，θ，z），將z軸作為對稱軸，此時(shí)應(yīng)力、應(yīng)變和位移都與θ方向無關(guān)，僅為r、z的函數(shù)，任一點(diǎn)的位移只有2個(gè)方向的分量，即沿r方向的徑向位移u和沿z方向的軸向位移w。

軸對稱問題的幾何關(guān)系為：

軸對稱問題的本構(gòu)關(guān)系為：

其中

由于采用多項(xiàng)式與徑向基函數(shù)相耦合的方法來處理本質(zhì)邊界條件，從而對于固體力學(xué)中的軸對稱問題，其平衡方程對應(yīng)的標(biāo)準(zhǔn)變分（弱）形式為：

其中，▽suT代表▽uT的對稱部分。將（13）式代入（14）式得控制方程為：

其中

（16）式中，

1.3 數(shù)值積分實(shí)現(xiàn)與求解流程

本文所有積分均是針對總體問題域Ω和總體表面力邊界而言的。為了計(jì)算這些總體積分，必須將該問題域離散成背景網(wǎng)格，背景網(wǎng)格高斯積分法雖然存在網(wǎng)格，但不參與插值過程，僅是為了完成區(qū)域積分和邊界積分，并不影響無網(wǎng)格法的實(shí)質(zhì)，從而不會削弱無網(wǎng)格法對于變動邊界或不連續(xù)面的處理，并且計(jì)算效率高，易于實(shí)現(xiàn)，因此本文為了得到（16）式的積分，采用獨(dú)立于場節(jié)點(diǎn)的背景網(wǎng)格，在背景網(wǎng)格的每一個(gè)積分子域（cell）內(nèi)采用高斯積分進(jìn)行計(jì)算。

依據(jù)以上思想，本文利用Fortran程序語言編程，實(shí)現(xiàn)了利用耦合多項(xiàng)式基的徑向點(diǎn)插值法進(jìn)行軸對稱問題的算法，計(jì)算流程如圖1所示。

圖1 計(jì)算流程圖

2 數(shù)值算例

一厚壁圓筒受均勻外壓示意如圖2所示，圓筒外壁半徑b＝2m，內(nèi)壁半徑為a＝1m，外壓p＝1.5kPa，E＝100GPa，υ＝0.3。

圖2 厚壁圓筒受均勻外壓

由于其對稱性，僅取1／4建模，對其進(jìn)行離散，離散后的節(jié)點(diǎn)分布如圖3所示，共布置節(jié)點(diǎn)253個(gè)，背景網(wǎng)格437個(gè)，在每個(gè)積分域內(nèi)采用3×3高斯積分。本文選擇MQ徑向基函數(shù)，其形參數(shù)取值為αc＝1.42，q＝1.63。

為檢驗(yàn)文中無網(wǎng)格法的計(jì)算效果，將計(jì)算結(jié)果和ANSYS有限元軟件的計(jì)算結(jié)果相比較，如圖4、圖5所示。在進(jìn)行ANSYS有限元軟件數(shù)值計(jì)算時(shí)，基于Plane42單元，節(jié)點(diǎn)數(shù)共計(jì)436個(gè)，單元數(shù)共計(jì)400個(gè)。

圖3 節(jié)點(diǎn)分布

圖4 沿壁厚方向的應(yīng)力計(jì)算結(jié)果與有限元解比較

圖5 徑向位移

為了更直觀地對比分析無網(wǎng)格法和有限元大型軟件ANSYS的數(shù)值計(jì)算結(jié)果，本文研究了厚壁圓筒受均勻外壓時(shí)1／4模型的應(yīng)力分布，其結(jié)果如圖6所示。

圖6 正應(yīng)力的計(jì)算結(jié)果與有限元解比較

從圖6可以看出，采用耦合多項(xiàng)式基的徑向點(diǎn)插值法計(jì)算所得應(yīng)力解與ANSYS有限元解吻合較好，取得了令人滿意的結(jié)果，同時(shí)也驗(yàn)證了該方法的有效性和合理性。

3 結(jié)束語

本文將耦合多項(xiàng)式基的徑向點(diǎn)插值無網(wǎng)格法推廣應(yīng)用于工程軸對稱問題中，并通過數(shù)值計(jì)算與有限元大型軟件ANSYS計(jì)算結(jié)果相比較，說明了耦合多項(xiàng)式基的徑向點(diǎn)插值法可以有效地求解工程軸對稱問題。數(shù)值結(jié)果表明，該無網(wǎng)格法計(jì)算精度高、穩(wěn)定性好，可直接施加本質(zhì)邊界條件，為將其應(yīng)用到軸對稱裂紋以及裂紋擴(kuò)展問題打下了基礎(chǔ)，具有較好的應(yīng)用前景。

本研究僅計(jì)算分析了厚壁圓筒受均勻擠壓的情況，對于擠壓凹模、受壓鋼瓶、柱形壓力容器、圓柱體鐓粗及多層包扎容器等工程上常見的軸對稱問題進(jìn)行簡化，建立無網(wǎng)格數(shù)值模型，也可以采用與此類似的方法進(jìn)行分析。該模型中形狀參數(shù)的選擇將會影響到該無網(wǎng)格法形函數(shù)的求解精度，為了獲得理想的結(jié)果，在進(jìn)行程序調(diào)試時(shí)，需對MQ函數(shù)中的參數(shù)做認(rèn)真仔細(xì)的選擇。此外，在相同規(guī)模下計(jì)算時(shí)，文中無網(wǎng)格法要比有限元法CPU耗時(shí)多一些。

［1］ 張 雄，劉 巖.無網(wǎng)格法［M］.北京：清華大學(xué)出版社，2004：1－6.

［2］ 趙光明.無單元法理論與應(yīng)用［M］.合肥：中國科學(xué)技術(shù)大學(xué)出版社，2010：1－14.

［3］ 錢德玲，程媛媛，邊燕飛.無網(wǎng)格法與有限元的耦合在線彈性斷裂力學(xué)中的應(yīng)用［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，32（7）：1061－1064.

［4］ 邊燕飛，何沛祥，程媛媛.一種自適應(yīng)影響域半徑無網(wǎng)格Galerkin法的應(yīng)用研究［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，32（4）：535－538.

［5］ Wang J G，Liu G R.A point interpolation meshless method based on radial basis functions［J］.J Numer Methods Eng，2002，54（11）：1623－1648.

［6］ 陳少林，李燕秀.一種高效的局部徑向基點(diǎn)插值無網(wǎng)格方法［J］.固體力學(xué)學(xué)報(bào)，2009，30（1）：100－105.

［7］ Liu G R，Gu Y T.A point interpolation method for two－dimensional solids［J］.International Journal for Numerical Methods in Engineering，2001，50（4）：937－951.

［8］ Wang J G，Liu G R.On the optimal shape parameters of radial basis functions used for 2－D meshless methods ［J］.Comput Methods Appl Mech Eng，2002，191：2611－2630.