基于sinh-sigmoidal組合式變換的薄型結(jié)構(gòu)邊界元分析方法*

張亞萍

(泰州職業(yè)技術(shù)學院 智能制造學院,江蘇 泰州 225300)

0 引 言

作為一種可拆卸的靜密封連接結(jié)構(gòu),墊片-法蘭-螺栓連接是石油管道及化工設(shè)備中必不可少的重要部件,目前已被廣泛應(yīng)用于機床裝備、航天航空等眾多工業(yè)領(lǐng)域[1-4]。

結(jié)構(gòu)的密封性能直接影響工業(yè)生產(chǎn)安全和企業(yè)生產(chǎn)利益。墊片作為連接結(jié)構(gòu)中的重要密封元件,其彈性力學性能也受到了普遍關(guān)注。

常見的密封方式是通過對上法蘭中的螺栓墊片施加預(yù)緊力。隨著所受預(yù)緊力的逐漸增大,墊片在受壓后,會表現(xiàn)出高度的非線性與復(fù)雜的時滯性。因此,準確評估初始預(yù)緊力及最大位移載荷下的壓縮回彈性能是保證其密封性能的重要前提[5]。

研究者們一直致力于厘清影響墊片-法蘭-螺栓連接結(jié)構(gòu)密封性能及連接可靠性的影響因素,并不斷優(yōu)化靜密封連接結(jié)構(gòu)設(shè)計的計算方法和相應(yīng)的試驗標準,以期建立起嚴謹有效的連接設(shè)計體系。歐美和日本等國家和地區(qū)已經(jīng)建立了基于泄漏率的墊片參數(shù)以及墊片壓縮回彈性能試驗方法,并形成了相應(yīng)的體系標準[6-8]。這套體系標準在測試墊片參數(shù)時,均假設(shè)施加在墊片上的應(yīng)力是均勻分布的。

然而,在實際的生產(chǎn)工作中,種種因素的影響會導(dǎo)致法蘭發(fā)生偏轉(zhuǎn),通過螺栓施加在墊片表面上沿徑向分布的接觸應(yīng)力是不均勻的。針對這一現(xiàn)象,現(xiàn)有文獻中普遍采用有限元法(FEM)對墊片的接觸應(yīng)力進行數(shù)值分析。

王璐等人[9]運用ABAQUS分析了當各螺栓施加的預(yù)緊力產(chǎn)生差異時,墊片接觸應(yīng)力的分布情況。盧志珍、朱樂等人[10-11]運用ANSYS自帶的襯墊單元模型Gasket,向墊片施加位移載荷,且步長為2 s,從而得出墊片的壓縮回彈性曲線,借此分析其密封性能。喻健良等人[12]運用ABAQUS,分析了不同加載方式和墊片型式對螺栓交互作用的影響。OMIYA Y等人[13-14]運用ABAQUS,分析了在各螺栓預(yù)緊力相同的情況下,墊片的接觸應(yīng)力分布情況,并通過實驗對此進行了驗證。

有限元法是目前一種主流的數(shù)值分析方法,并已經(jīng)得到了廣泛的應(yīng)用。但在計算薄壁結(jié)構(gòu)、應(yīng)力集中、無限域等問題時,其計算精度及計算效率較低。為彌補有限元法的不足,近年來諸如無網(wǎng)格法[15]、邊界元法[16-19]等數(shù)值分析方法取得了長足進步,各個瓶頸問題也相繼得到了解決。

邊界元法(BEM)是在邊界積分方程理論基礎(chǔ)上建立起來的一種純邊界離散型數(shù)值分析方法。相較于其他數(shù)值分析方法,其具有離散單元數(shù)量較少、數(shù)據(jù)準備簡單、計算量小的優(yōu)勢。但利用BEM分析墊片等薄型結(jié)構(gòu)的彈性力學問題時,求解的大量近奇異積分會出現(xiàn)近奇異性,且被積函數(shù)會隨距離的減小而急劇變化,需要大量高斯積分點才能獲得較好的計算精度。這一情況嚴重影響了該方法的計算效率。

針對薄型結(jié)構(gòu)分析中的積分近奇異性消除問題,SLADEK V等人[20]于1993年提出了正則化法。2014年,ZHANG Xi-yu等人[21]也提出了一種基于薄型結(jié)構(gòu)混合邊界節(jié)點法的理論分析方法,數(shù)值算例表明了該方法具有一定的適用性。2019年,張見明等人[22]提出了一種新的單元細分法(基于體二叉樹),大幅提高了奇異和近奇異積分的計算穩(wěn)定性及精度。

在不同的方法中,非線性變換法能更好地利用計算機軟件,對其進行模塊化編程,并且可以不依賴積分核函數(shù);但其計算精度和穩(wěn)定性受到積分域內(nèi)部投影點位置的影響。

筆者將針對弱化投影點影響力的新型非線性變換法進行研究,提出一種非線性組合式sinh-Sigmoidal變換方法,并將其集成到邊界元程序,以消除積分在徑向和角度方向的近奇異性;然后,分別在理想工況和實際工況下,采用BEM和FEM程序分析薄型墊片的彈性靜力學特性,以驗證該方法的有效性。

1 邊界元法

筆者以彈性力學問題為例進行說明,將邊界元法應(yīng)用于實際力學問題的求解。

經(jīng)典彈性力學問題的平衡方程和邊界條件描述如下[23]:

(1)

其中:對于三維問題,均采用指標形式記法,下標i,j=1,2,3表示對其所代表的變量求偏導(dǎo)數(shù),即σij,j=?σij/?xj。

同一變量中,若下標重復(fù)出現(xiàn)即表示求和,其幾何方程和物理方程如下:

(2)

式中:δij為Kronecker符號,當下標i=j時通常取值為1,否則為零;λ,μ均為Lame彈性常數(shù)。

λ,μ與彈性模量E、泊松比ν以及剪切模量G之間的關(guān)系如下:

(3)

經(jīng)典彈性力學問題相應(yīng)的邊界積分方程常表示如下:

(4)

對于三維各向同性體,其基本解為開爾文(Kelvin)解。位移和面力的基本解表示如下[24]:

(5)

(6)

式中:r(x,y)為源點y與場點x間的歐幾里得距離。

其中:位移的基本解滿足愛因斯坦求和約定,當l=i時,δli=1,其余情況取值均為零。

假設(shè)源點y位于光滑的積分邊界上,則取系數(shù)clj(y)=1/2。在不考慮體力影響的情況下,可得只包含邊界積分項的積分方程,通過單元插值法將所求解的具體問題的邊界進行離散化處理。

對位移和面力進行等參差值計算如下:

(7)

(8)

式中:m為邊界離散的單元總數(shù)量;Nk(ξ,η)則為離散單元內(nèi)部第k個插值點所對應(yīng)的形函數(shù)。

筆者通過廣義達菲變換,將坐標轉(zhuǎn)換至(ξ,η)空間,得到了正則化的參數(shù)坐標。在給定的相應(yīng)邊界條件下,將場點x取遍積分單元內(nèi)所有插值點的位置,可得3N×3N規(guī)模的線性方程組(其中,N為總的插值點數(shù)量)。

通過計算機編程求解,筆者可得邊界上任意未知的位移或面力值。而彈性體內(nèi)任意一點的位移、應(yīng)力、應(yīng)變的值則可以根據(jù)式(4)計算得出。

2 sinh-sigmoidal組合式變換法

筆者利用BEM分析薄型結(jié)構(gòu)的彈性力學問題,求解大量的近奇異積分,公式如下:

(9)

式中:f為不存在近奇異性的光滑函數(shù);x為域內(nèi)被積點;y為域外源點;r為兩點間的距離。

當分母rl趨近于零時,積分在徑向和角度方向均會出現(xiàn)近奇異性,被積函數(shù)隨距離的減小而急劇變化,需要大量高斯積分點才能獲得較好的計算精度,但其嚴重影響了計算效率。

因此,精確高效地求解近奇異積分是利用BEM分析薄型結(jié)構(gòu)問題的關(guān)鍵?；谧赃m應(yīng)分塊技術(shù)和不同坐標變換,筆者提出一種非線性sinh-sigmoidal組合變換法,將其與原始極坐標變換下(ρ,θ)坐標空間中的近奇異積分相結(jié)合,可以在提高積分計算精度的同時大幅降低計算量,彌補傳統(tǒng)變換法的不足。

2.1 原始極坐標變換

局部坐標系中的距離函數(shù)如圖1所示。

圖1 局部坐標系中的距離函數(shù)

在三角形單元中,距離源點y最近的點為投影點xc,兩者的連線垂直于積分單元中xc處的切線。

根據(jù)源點y的位置,即可確定投影點xc的位置。在投影點附近使用泰勒公式展開如下[25-26]:

(10)

式中:(ξ0,η0)為投影點在局部空間坐標系(ξ,η)內(nèi)的參數(shù)坐標。

以(ξ0,η0)為極坐標原點變換所得,則ξ、η的表達式如下:

(11)

將式(11)代入式(10)可得:

(12)

從源點y到場點x的實際距離r又可以表達如下:

(13)

其中:ω(θ)=r0/A(θ),A(θ)=|Ak(θ)|。

在局部坐標空間系(ξ,η)中,筆者對近奇異積分的一般形式即式(9)進行原始極坐標變換,其表達式如下:

(14)

由式(14)可得,原始極坐標變換將積分空間從(ξ,η)轉(zhuǎn)換至(ρ,θ),且徑向和角度方向上均出現(xiàn)了積分近奇異性。

2.2 徑向變換

為保證投影點位于積分形狀良好的子三角形單元內(nèi)部,降低積分形狀對計算精度的影響。筆者對每個積分單元運用自適應(yīng)分塊技術(shù),將sinh變換重新定義如下:

r(s)=ω(θ)sinh(μ1s-η1)

(15)

將積分區(qū)間[0,ρ(θ)]和[θm,θm+1]分別轉(zhuǎn)換至區(qū)間[-1,1],可得:

(16)

因為投影點位于積分單元邊界的高階近奇異積分精度不夠,所以迭代sinh變換是必要的。其公式如下:

s(u)=a+bsinh(μ2u-η2)

(17)

其中:a=-1,b=π/2μ1。將引入的新變量的積分區(qū)間轉(zhuǎn)換至具體區(qū)間[-1,1],簡化編程過程,可得下式:

(18)

將式(15)和式(17)代入式(14),可得:

(19)

筆者將上述積分區(qū)間從未知量轉(zhuǎn)換成已知量,可消除積分核函數(shù)在徑向上的積分近奇異性;但角度方向的積分近奇異性依然存在。

2.3 角度方向近奇異性

在典型的sigmoidal變換中,其核函數(shù)能有效地向積分區(qū)間兩端大量聚集高斯積分點。利用這個特點,可將其應(yīng)用于消除角度方向積分近奇異性。常用變換式如下:

(20)

式中:w為變換變量的冪指數(shù),且w≥1。

由此又可得到:

(21)

為了消除積分的雙向近奇異性,并提升計算的精度,將式(21)代入式(19),可得:

(22)

在上述轉(zhuǎn)換過程中,通常取w=2,以獲得滿意的精度。

3 薄型墊片彈性靜力學特性分析

3.1 理想工況邊界元分析

根據(jù)上下表面及內(nèi)外環(huán)設(shè)定的邊界條件,以及無體力項的彈性力學邊界積分方程,筆者對墊片所有邊界節(jié)點處面力的變化情況進行多方向分析,并通過不斷加大網(wǎng)格劃分密度和增加計算節(jié)點數(shù)量,以此來驗證邊界元法的收斂性。

筆者對該薄型墊片的各個面施加x、y、z方向的位移邊界約束條件,如下:

(23)

為了驗證邊界元法計算所得結(jié)果的收斂性,筆者取4種不同的網(wǎng)格尺寸以及單元、節(jié)點數(shù)量,對薄型墊片所有邊界計算節(jié)點進行x、y、z方向上面力變化情況的誤差分析。

邊界節(jié)點各方向面力的平均誤差計算結(jié)果如表1所示。

表1 邊界節(jié)點各方向面力的平均誤差計算

邊界節(jié)點各方向面力的平均誤差值分析結(jié)果如圖2所示。

圖2 邊界節(jié)點各方向面力的平均誤差值分析

圖2中結(jié)果顯示:隨著網(wǎng)格尺寸的不斷細分,單元和節(jié)點數(shù)量的不斷增加,計算節(jié)點處的誤差值不斷減小。

運用BEM程序計算可得,x、y、z方向上節(jié)點位移的數(shù)值結(jié)果能有效地收斂到精確解,說明該方法在分析薄壁結(jié)構(gòu)彈性力學問題時是精確可行的。

3.2 實際工況墊片受力變形分析

上述結(jié)果已經(jīng)驗證了邊界元法在理想工況下(在各方向邊界條件均可以用公式描述),分析薄型墊片所得結(jié)果的有效性與正確性。

為進一步驗證該方法在實際工況下同樣具有可行性與一定的精度,筆者分別采用兩種方法(BEM、FEM)分析計算墊片受外力影響下的變形情況。

薄型墊片的約束條件及載荷施加如圖3所示。

圖3 墊片模型受力情況

因薄型墊片的結(jié)構(gòu)輕量化,此處不考慮自重產(chǎn)生的變形。筆者將其一端面固定,另一端面施加壓力。

邊界元法程序網(wǎng)格劃分如圖4所示。

圖4 邊界元法網(wǎng)格劃分

為了獲得更為精確的變形結(jié)果,筆者在邊界元法分析中使用了1 296個三角形單元,域內(nèi)及表面共計使用2 113個節(jié)點。

筆者采用ABAQUS有限元軟件對墊片模型進行網(wǎng)格劃分,其結(jié)果如圖5所示。

圖5 有限元法網(wǎng)格劃分

在有限元法分析中,筆者共計使用了126 350個線性六面體單元(C3D8R)以及142 791個節(jié)點。

采用ABAQUS有限元軟件對模型進行分析,可以得到z方向的位移變形云圖,如圖6所示。

圖6 基于有限元法的z方向位移云圖

墊片在z方向位移變形范圍區(qū)間為[-2.069,0]。

為進一步驗證邊界元法分析的精確性和高效性,筆者在其厚度方向均勻選取10個取樣點,并進行對比分析。

樣點坐標及計算結(jié)果對比數(shù)據(jù)如表2所示。

表2 樣點坐標及2種方法計算結(jié)果對比

由表2數(shù)據(jù)可知:采用邊界元法,以極少的單元和節(jié)點數(shù)量計算所得z方向位移的變形結(jié)果,與采用有限元法分析所得到的結(jié)果相近。

厚度方向FEM與BEM的位移計算結(jié)果對比如圖7所示。

圖7 厚度方向FEM與BEM的位移計算結(jié)果對比

由圖7可知:運用BEM程序計算得到的樣點處z方向位移線性圖,與ABAQUS軟件計算所得圖形趨勢基本吻合,并且從10個取樣點的具體數(shù)據(jù)值而言,上述2種方法間的計算誤差均小于1%。

由此可見,該算例也充分驗證了邊界元法在實際工況下分析薄型墊片彈性力學問題的高效性以及精確性。

4 結(jié)束語

針對薄型結(jié)構(gòu)邊界元數(shù)值分析過程中的近奇異積分計算困難問題,筆者提出了一種基于sinh-Sigmoidal組合變換的解決方案,消除了積分在徑向和角度方向的近奇異性,并提高了邊界元法的計算效率和計算精度。

通過將其與有限元分析方法(FEM)進行比較,結(jié)果驗證了該方法的有效性,其是一種適用于薄型墊片彈性靜力學問題的有效數(shù)值分析方法。

研究結(jié)論如下:

1)采用迭代sinh變換可消除積分的徑向奇異性,采用Sigmoidal變換可消除角度方向奇異性;

2)通過理想工況下的薄型墊片算例分析,驗證了BEM在計算墊片結(jié)構(gòu)面力時的收斂性和正確性;

3)通過實際工況下的薄型算例分析,對比了10個取樣點處BEM和FEM的位移計算結(jié)果,計算誤差均小于1%;但BEM方法的計算效率優(yōu)勢明顯,驗證了BEM計算的高效性以及精確性;

4)在薄型墊片的設(shè)計中,可利用BEM對墊片進行彈性力學分析,得到在表面預(yù)緊力作用下的變形曲線,再進行相應(yīng)的補償設(shè)計,以保證薄型墊片的連接可靠性,提高其密封性能。

當前的研究主要考慮的是投影點位于積分單元內(nèi)部的情況。對投影點位于外部的情況,有待筆者進一步深入研究。同時,由于在當前研究中,所分析對象為結(jié)構(gòu)簡單的圓形薄墊片,要使BEM方法能得到廣泛的應(yīng)用,未來還需針對更復(fù)雜的結(jié)構(gòu)開展研究。