一種與時(shí)間步長(zhǎng)相關(guān)的奇異單元細(xì)分法

李 源， 張見(jiàn)明， 鐘玉東， 千紅濤

(1.河南師范大學(xué) 計(jì)算機(jī)與信息工程學(xué)院，河南 新鄉(xiāng) 453007； 2.湖南大學(xué) 機(jī)械與運(yùn)載工程學(xué)院，湖南 長(zhǎng)沙 410082； 3.河南工學(xué)院 機(jī)械工程系，河南 新鄉(xiāng) 453003)

0 引言

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，數(shù)值分析與仿真技術(shù)已經(jīng)被廣泛應(yīng)用于工業(yè)設(shè)計(jì)和生產(chǎn)過(guò)程中，例如靜力學(xué)分析[1]、熱力學(xué)分析[2]、裂紋擴(kuò)展分析[3]等，數(shù)值計(jì)算方法是有效的分析手段之一.邊界元法[4]利用微分方程的解析基本解作為邊界積分方程的核函數(shù)，是一種半解析的數(shù)值計(jì)算方法，通常精度較高，特別是對(duì)于求解裂紋等奇異性問(wèn)題[5]有天然優(yōu)勢(shì).但奇異積分則是邊界元法求解物理問(wèn)題時(shí)的難點(diǎn)，當(dāng)積分場(chǎng)點(diǎn)和源點(diǎn)位于同一個(gè)單元內(nèi)時(shí)，兩者距離可能為零，由于距離因子r位于核函數(shù)的分母上，因此r趨向于大于0時(shí)將導(dǎo)致核函數(shù)值無(wú)窮大，得到錯(cuò)誤的結(jié)果，因此需要采取措施避免場(chǎng)點(diǎn)和源點(diǎn)的重合.在傳統(tǒng)的單元細(xì)分法中，根據(jù)奇異單元的形狀以及源點(diǎn)在單元中的位置，將源點(diǎn)與單元中各節(jié)點(diǎn)相連接，原單元被細(xì)分為若干個(gè)三角形子單元，從而能夠避免源點(diǎn)和場(chǎng)點(diǎn)的重合[6].Zhang等[7-8]曾將該方法與一種新的坐標(biāo)變換法相結(jié)合消除了三維勢(shì)問(wèn)題的奇異性，后來(lái)又提出了一種自適應(yīng)單元細(xì)分法[9]來(lái)提高積分精度.周煥林等[10]提出了一種解析積分算法和單元細(xì)分相結(jié)合的方法解決了二維各向異性位勢(shì)問(wèn)題中的近奇異積分.

上述細(xì)分法應(yīng)用于靜態(tài)問(wèn)題時(shí)計(jì)算精度較高. 對(duì)于動(dòng)態(tài)問(wèn)題，例如筆者所關(guān)注的彈性動(dòng)力學(xué)問(wèn)題[11-12]基本解更加復(fù)雜，不僅包含空間變量，還包含時(shí)間變量，且為分段函數(shù)(表現(xiàn)為縱波和橫波兩個(gè)波動(dòng)前沿).如果仍然采用傳統(tǒng)的單元細(xì)分法，則無(wú)法將波動(dòng)前沿對(duì)核函數(shù)的影響考慮進(jìn)去，奇異積分的精度必然會(huì)受到影響，因此筆者根據(jù)三維問(wèn)題中彈性波傳播的規(guī)律以及核函數(shù)的分段特性，將波速和時(shí)間步長(zhǎng)考慮在內(nèi)，提出了一種新的單元細(xì)分法.結(jié)合文獻(xiàn)[8]中的坐標(biāo)變換有效地消除了基本解中的奇異性，提高了積分的精度，最后通過(guò)兩個(gè)算例驗(yàn)證該方法的準(zhǔn)確性和計(jì)算精度.

1 彈性動(dòng)力學(xué)問(wèn)題的時(shí)域邊界積分方程

對(duì)于均勻、各向同性的線彈性材料而言，彈性動(dòng)力學(xué)問(wèn)題的位移運(yùn)動(dòng)方程(亦稱為控制方程)如下[12]：

(1)

(2)

式中：E表示彈性模量；ν表示泊松比.

由于筆者研究的重點(diǎn)為奇異積分處理中的單元細(xì)分方法，而奇異積分主要出現(xiàn)在邊界積分中，因此為了敘述方便，設(shè)問(wèn)題的初始條件和體積力均為零(不為零時(shí)關(guān)于奇異積分的處理仍然是一樣的)，此時(shí)控制方程(1)對(duì)應(yīng)的時(shí)域邊界積分方程為[12]：

(3)

2 邊界積分方程中奇異積分的處理

方程(3)中的基本解經(jīng)過(guò)時(shí)間離散和積分后稱為核函數(shù)，根據(jù)該問(wèn)題位移核函數(shù)的表達(dá)式作出能表示其變化趨勢(shì)的示意圖,如圖1所示.圖1中在r=c1Δt和r=c2Δt處均發(fā)生了函數(shù)值的突變(圖1中：橫坐標(biāo)r代表場(chǎng)點(diǎn)距離源點(diǎn)的距離；縱坐標(biāo)Uij代表位移核函數(shù)值；c1代表縱波的波速；c2代表橫波的波速；Δt代表時(shí)間步長(zhǎng))，即在波動(dòng)前沿處核函數(shù)值會(huì)產(chǎn)生突變.當(dāng)源點(diǎn)位于積分單元內(nèi)且c1Δt

圖1 位移核函數(shù)曲線示意圖Fig.1 Diagram of the displacement kernel function curve

圖2 波動(dòng)前沿示意圖Fig.2 Diagram of the wave front edge

基于以上考慮，根據(jù)單元尺寸和時(shí)間步長(zhǎng)之間的關(guān)系，筆者將奇異積分的單元分3種情況進(jìn)行細(xì)分.

(1)當(dāng)c2Δt≥d時(shí)，細(xì)分方法如圖3所示，傳統(tǒng)單元細(xì)分法采用的就是這種方式，而在本方法中該細(xì)分方式只是其中之一.當(dāng)源點(diǎn)位于單元的角點(diǎn)上時(shí)，將源點(diǎn)與對(duì)角點(diǎn)相連，如圖3(a)所示，原單元被細(xì)分為兩個(gè)三角形子單元；當(dāng)源點(diǎn)位于單元的邊上時(shí)，如圖3(b)所示，將源點(diǎn)與其他兩個(gè)節(jié)點(diǎn)分別相連形成3個(gè)三角形子單元；當(dāng)源點(diǎn)位于單元內(nèi)時(shí)，如圖3(c)所示，將源點(diǎn)與4個(gè)節(jié)點(diǎn)分別相連形成4個(gè)三角形子單元.在細(xì)分后的子單元內(nèi)積分時(shí)，不會(huì)再出現(xiàn)場(chǎng)點(diǎn)與源點(diǎn)重合的情況.第1種細(xì)分法出現(xiàn)的概率很小，因?yàn)闀r(shí)間步長(zhǎng)過(guò)大時(shí)計(jì)算結(jié)果的振幅衰減過(guò)快，精度很低.

圖3 第1種單元細(xì)分方式Fig.3 The first kind of element subdivision method

(2)當(dāng)c2Δt

圖4 第2種單元細(xì)分方式Fig.4 The second kind of element subdivision method

(3)當(dāng)d≥c1Δt時(shí)，細(xì)分方法如圖5所示.此時(shí)橫波和縱波的波動(dòng)前沿都與單元相交，如圖中圓弧形點(diǎn)畫(huà)線所示.首先按照第二種情況的步驟處理包含橫波波動(dòng)前沿的四邊形塊；然后求出縱波波動(dòng)前沿與單元邊界的交點(diǎn)，根據(jù)交點(diǎn)位置將單元的其他部分劃分為若干個(gè)四邊形塊，最終得到的分塊結(jié)果如圖5中的實(shí)線所示.同樣的，只在包含源點(diǎn)的三角形塊中采用奇異積分即可，其他分塊區(qū)域只需采用近奇異積分或正則積分.

圖5 第3種單元細(xì)分方式Fig.5 The third kind of element subdivision method

與傳統(tǒng)的單元細(xì)分法相比，這種改進(jìn)的方法能夠根據(jù)波動(dòng)前沿的位置動(dòng)態(tài)改變單元細(xì)分方式，在被波動(dòng)前沿波及的區(qū)域布置更多的積分點(diǎn)，得到更高的奇異積分精度，而在波動(dòng)前沿還未到達(dá)的區(qū)域無(wú)需布置積分點(diǎn).需要注意的是，奇異積分只發(fā)生在第一個(gè)分析步，后面的分析步無(wú)需使用單元細(xì)分法.下面將通過(guò)兩個(gè)數(shù)值算例對(duì)比改進(jìn)后和改進(jìn)前計(jì)算結(jié)果的精度.

3 數(shù)值算例

3.1 懸臂梁在Heaviside類型載荷下的動(dòng)態(tài)響應(yīng)

本例以懸臂梁為分析模型，模型的幾何尺寸如圖6所示.邊界條件：左端面約束所有方向的自由度，右端面施加值為1 000 Pa的均布載荷，載荷為沿x軸正方向的Heaviside類型的均布力.材料參數(shù)：彈性模量E=1.1×105Pa，密度ρ=2.0 kg/m3，泊松比ν=0.當(dāng)泊松比為0時(shí)，該問(wèn)題可通過(guò)理論方法計(jì)算出解析解.因此本例中將泊松比設(shè)置為0，是為了得到數(shù)值解的計(jì)算誤差，在程序中仍然將該問(wèn)題當(dāng)作三維問(wèn)題來(lái)計(jì)算.根據(jù)文獻(xiàn)[13]可得到懸臂梁上任意一點(diǎn)任意時(shí)刻的位移響應(yīng)計(jì)算公式：

(4)

圖6 懸臂梁模型和載荷曲線示意圖Fig.6 Diagram of the model and load curve of the cantilever beam

采用尺寸為1 m的線性四邊形單元離散模型.首先取時(shí)間步長(zhǎng)Δt=0.003 56 s，得到前10個(gè)分析步的計(jì)算結(jié)果.表1給出了懸臂梁自由端位移的響應(yīng)，表中最后一列的數(shù)值等于改進(jìn)前結(jié)果的相對(duì)誤差減去改進(jìn)后結(jié)果的相對(duì)誤差，可以看出改進(jìn)后相對(duì)誤差明顯減小，尤其是存在奇異積分的第一個(gè)分析步結(jié)果的相對(duì)誤差減小了15.5%.將結(jié)果繪制成曲線，如圖7所示.圖7中exact曲線代表精確解，traditional曲線代表改進(jìn)前的結(jié)果，improved曲線代表改進(jìn)后的結(jié)果.可以明顯看出改進(jìn)后的結(jié)果與精確解更加吻合.

表1 Δt=0.003 56 s時(shí)懸臂梁自由端位移隨時(shí)間的變化Tab.1 The displacement of the free end of the cantilever beam with time when Δt=0.003 56 s

圖7 Δt=0.003 56 s時(shí)改進(jìn)前后結(jié)果的對(duì)比曲線Fig.7 The comparison curves of the results before and after the improvement when Δt=0.003 56 s

將時(shí)間步長(zhǎng)減小為原來(lái)的1/2,Δt=0.001 78 s，圖8給出了該步長(zhǎng)下改進(jìn)前后計(jì)算結(jié)果的對(duì)比.從圖中可以看出，改進(jìn)前的結(jié)果在某些時(shí)刻會(huì)出現(xiàn)不穩(wěn)定的現(xiàn)象，改進(jìn)后的結(jié)果與精確解更為吻合，再次驗(yàn)證了本文方法的有效性.

圖8 Δt=0.001 78 s時(shí)改進(jìn)前后結(jié)果的對(duì)比曲線Fig.8 The comparison curves of the results before and after the improvement when Δt=0.001 78 s

3.2 帶孔平板在動(dòng)態(tài)載荷下的響應(yīng)

該問(wèn)題的分析模型如圖9(a)所示，圖9(a)中幾何尺寸單位為mm，邊界條件：下端面約束所有方向的自由度，上端面施加隨時(shí)間變化的均布載荷，載荷曲線如圖9(b)所示，初始時(shí)刻載荷為零，在t=3.75×10-5s時(shí)達(dá)到最大值118.125 MPa.模型的材料參數(shù)為：彈性模量E=6.9×104MPa，密度ρ=2.7×10-9kg/mm3，泊松比ν=0.3.

圖9 帶孔平板的幾何模型及載荷曲線示意圖Fig.9 The geometric model and load curve of a plate with a hole

圖10 t=2.5×105 s時(shí)帶孔平板的位移分布云圖Fig.10 The displacement nephogram of the plate with a hole when t=2.5×105 s

采用尺寸為10 mm的線性三角形單元離散模型，該算例無(wú)解析解，將本文的計(jì)算結(jié)果與有限元分析軟件ABAQUS的結(jié)果進(jìn)行對(duì)比，為了得到更加準(zhǔn)確的對(duì)比數(shù)據(jù)，將有限元網(wǎng)格尺寸取為5 mm.

取t=2.5×10-5s時(shí)所有節(jié)點(diǎn)的位移值繪制成位移分布云圖，并與有限元位移云圖進(jìn)行比較，如圖10所示，本文方法的結(jié)果云圖和有限元法基本一致，再次驗(yàn)證了該方法的正確性.

為了觀察小孔附近的應(yīng)力集中現(xiàn)象，根據(jù)t=2.5×10-5s時(shí)所有節(jié)點(diǎn)的mises應(yīng)力值繪制成應(yīng)力分布云圖，并與有限元結(jié)果云圖比較，如圖11所示，應(yīng)力分布與有限元的分析結(jié)果基本一致，可以明顯看出小孔附近的應(yīng)力值最大.

圖11 t=2.5×10-5 s時(shí)帶孔平板的應(yīng)力分布云圖Fig.11 The stress nephogram of the plate with a hole when t=2.5×10-5 s

4 結(jié)論

筆者針對(duì)瞬態(tài)彈性動(dòng)力學(xué)問(wèn)題，采用時(shí)域邊界元法進(jìn)行數(shù)值求解，給出了該問(wèn)題時(shí)域邊界積分方程的一般形式.通過(guò)研究基本解函數(shù)的性質(zhì)，提出了一種適用于動(dòng)態(tài)問(wèn)題的與時(shí)間步長(zhǎng)相關(guān)的奇異單元細(xì)分法.在本研究中，首先由時(shí)間步長(zhǎng)和波速計(jì)算出兩個(gè)波動(dòng)前沿的傳播半徑，然后結(jié)合源點(diǎn)的位置將奇異積分單元分3種情況進(jìn)行細(xì)分，細(xì)分后再通過(guò)坐標(biāo)變換法消除弱奇異性，剛體位移法消除強(qiáng)奇異性.與只考慮源點(diǎn)位置的傳統(tǒng)細(xì)分法相比，該方法提高了奇異積分的精度，從而也提高了最終計(jì)算結(jié)果的精度.最后的兩個(gè)算例也驗(yàn)證了方法的有效性，結(jié)果顯示：對(duì)于存在奇異積分的第一個(gè)分析步，結(jié)果誤差能夠減小15.5%；利用該方法得到的仿真結(jié)果能夠準(zhǔn)確模擬位移和應(yīng)力的分布與變化.