機(jī)械系統(tǒng)靜態(tài)特性的有限元耦合分析研究

劉 潔

（延安大學(xué)化學(xué)與化工學(xué)院，陜西 延安716000）

0 引言

商用有限元軟件是一個(gè)黑盒子，其二次開發(fā)只能做一些外圍工作，其核心接觸算法屬于商業(yè)秘密。現(xiàn)在需要修改接觸算法及其收斂條件，才能引入結(jié)合面的特性參數(shù)進(jìn)行耦合分析，并且需要設(shè)計(jì)新的算法才能求解整機(jī)。基于有限元法，研究重點(diǎn)是把接觸算法和結(jié)合部理論進(jìn)行融合，設(shè)計(jì)耦合求解算法，并編制耦合分析軟件，最后用實(shí)驗(yàn)進(jìn)行驗(yàn)證。

1 有限元耦合分析

1.1 建立接觸系統(tǒng)的求解方程

設(shè)機(jī)械系統(tǒng)Ⅰ與機(jī)械系統(tǒng)Ⅱ相互接觸，每個(gè)機(jī)械系統(tǒng)可以是單個(gè)零件，也可以是由許多零部件通過柔性結(jié)合部連接成的裝配體［1］。建立機(jī)械系統(tǒng)Ⅰ與Ⅱ的離散平衡方程式（1）和（2），并對(duì)其系數(shù)剛度矩陣進(jìn)行分塊，然后引入已知位移和力邊界條件。

U1為系統(tǒng)Ⅰ非接觸邊界節(jié)點(diǎn)的位移列陣；U2為系統(tǒng)Ⅰ接觸邊界節(jié)點(diǎn)的位移列陣；F1為系統(tǒng)Ⅰ非接觸邊界節(jié)點(diǎn)的力邊界條件列陣；F2為系統(tǒng)Ⅰ接觸邊界節(jié)點(diǎn)的接觸力列陣。

U3為系統(tǒng)Ⅱ接觸邊界節(jié)點(diǎn)的位移列陣；U4為系統(tǒng)Ⅱ非接觸邊界節(jié)點(diǎn)的位移列陣；F3為系統(tǒng)Ⅱ接觸邊界節(jié)點(diǎn)的力邊界條件列陣；F4為系統(tǒng)Ⅱ非接觸邊界節(jié)點(diǎn)的接觸力列陣。

接觸力滿足F2＝－F3，但是接觸力的大小是未知的，所以方程無法直接求解。

1.2 結(jié)合部子結(jié)構(gòu)的罰函數(shù)法有限元

將接觸區(qū)域的非嵌入條件以及其他附加條件作為懲罰項(xiàng)引入接觸域子結(jié)構(gòu)的泛函中，將接觸問題轉(zhuǎn)化為關(guān)于接觸剛度的優(yōu)化問題［2－4］。接觸剛度的物理意義是結(jié)合面抵抗變形的能力，相當(dāng)于在接觸點(diǎn)對(duì)之間嵌入1個(gè)法向、2個(gè)切向壓縮彈簧，來抵抗節(jié)點(diǎn)間法向和切向的相對(duì)位移。

由式（1）和式（2）經(jīng)過靜力凝聚消去U1，U4后得到：

建立結(jié)合部子結(jié)構(gòu)的泛函：

∏U為接觸域子結(jié)構(gòu)不包括接觸條件的總勢(shì)能；∏CP為用罰函數(shù)法引入接觸約束條件的附加泛函；KJ為接觸剛度；JJ2，JJ3為總體坐標(biāo)系與結(jié)合部局部坐標(biāo)系之間的坐標(biāo)變換矩陣。

根據(jù)最小勢(shì)能原理，在滿足邊界條件的可能位移中，真實(shí)位移滿足非嵌入條件，并使系統(tǒng)總勢(shì)能取得極小值［5］。也就是接觸的最終狀態(tài)由系統(tǒng)平衡方程、力和位移邊界條件以及接觸域邊界的非嵌入條件所惟一確定。

當(dāng)泛函∏C取極值時(shí)，即

得到：

整理合并后得到接觸域子結(jié)構(gòu)的控制方程為：

由罰函數(shù)法得到的顯式方程式（6），可以采用基于載荷增量迭代的中心差分法求解。選用足夠大等間距的罰因子數(shù)列δ1，δ2，δ3，…，δn，對(duì)于某個(gè)給定的δk，當(dāng)接觸表面層的最大相對(duì)位移量λ2max≤ε（ε為給定的最大相對(duì)位移量）時(shí)，認(rèn)為迭代收斂。

罰函數(shù)法的優(yōu)點(diǎn)是算法簡單，易于編程［6］。但接觸條件只是近似被滿足，在接觸域存在穿透，最大相對(duì)位移量由給定的常量ε來控制。當(dāng)選用的罰因子δk過大時(shí)，此時(shí)的接觸剛度會(huì)很大，從而引起計(jì)算結(jié)果的病態(tài)。

1.3 對(duì)接觸域的精確處理

如果對(duì)計(jì)算精度的要求比較高時(shí)，需要檢驗(yàn)罰函數(shù)法的計(jì)算結(jié)果是否合理，并做進(jìn)一步精確處理。

設(shè)第k個(gè)接觸點(diǎn)對(duì)之間由接觸表面層變形產(chǎn)生的法向位移為λ2k，切向相對(duì)位移為λ1k，λ3k，則

合并寫成列陣得到：

則第k個(gè)接觸點(diǎn)對(duì)之間的接觸面壓為：

將非線性接觸面壓等效為節(jié)點(diǎn)載荷得：

［N］為單元形函數(shù)矩陣。

將結(jié)合部的相對(duì)變形｛λk｝和接觸力｛FC｝進(jìn)行集成裝配，得到總的列陣｛λ｝和｛FC｝。

聯(lián)立式（6）和式（12），使用增量迭代法求解，最終讓接觸表面層的最大相對(duì)位移量等于接觸表面層的相對(duì)變形。當(dāng)滿足收斂條件后，由式（12）得到接觸節(jié)點(diǎn)的接觸力FC列陣，將FC帶入式（1）和式（2）即可求解機(jī)械系統(tǒng)Ⅰ與Ⅱ。

一般情況下，機(jī)械系統(tǒng)Ⅰ沒有足夠的位移邊界條件，其剛度矩陣具有奇異性。還需要在機(jī)械系統(tǒng)Ⅰ上引入一個(gè)位移約束，并推出此時(shí)的剛體位移公式和由此產(chǎn)生的整體附加平衡條件，與式（1）聯(lián)立求解。

機(jī)械系統(tǒng)Ⅰ和Ⅱ都應(yīng)滿足整體附加平衡條件：

因結(jié)合部相對(duì)位移的數(shù)量級(jí)一般小于等于微米級(jí)，而接觸區(qū)域的零件宏觀變形和本體結(jié)構(gòu)及外載荷相關(guān)，進(jìn)行耦合時(shí)有以下2種迭代方案。

a.實(shí)時(shí)處理。當(dāng)接觸區(qū)域的零件宏觀變形和結(jié)合部的相對(duì)變形處于同一個(gè)數(shù)量級(jí)（或者相差較小）時(shí)，則在每一次迭代過程中均代入結(jié)合面的基礎(chǔ)特性參數(shù)，對(duì)接觸域進(jìn)行實(shí)時(shí)修正。

b.最終處理。當(dāng)接觸區(qū)域的零件宏觀變形和結(jié)合部的相對(duì)變形相差幾個(gè)數(shù)量級(jí)時(shí)，先不代入結(jié)合面的基礎(chǔ)特性參數(shù)，即先用罰函數(shù)法求出大致的接觸域后，再代入結(jié)合面的基礎(chǔ)特性參數(shù)做準(zhǔn)確處理，對(duì)接觸域進(jìn)行精確的修正。

2 誤差分析

誤差分析也是有限元計(jì)算中的重要問題之一，在計(jì)算過程中需要估計(jì)誤差，如果誤差過大，則需要改用高階單元，或者增加局部的單元密度，來保證計(jì)算精度。

2.1 計(jì)算誤差

2.1.1 舍入誤差

由于計(jì)算機(jī)字長所限，需要對(duì)原始數(shù)據(jù)、中間結(jié)果和最終結(jié)果保留有限位數(shù)字，有限位小數(shù)后的尾數(shù)部分做四舍五入處理，這種由四舍五入產(chǎn)生的誤差被稱為舍入誤差。舍入誤差是不可避免的，而且對(duì)整體計(jì)算精度的影響非常小。

2.1.2 截?cái)嗾`差

由于實(shí)際運(yùn)算只能完成有限項(xiàng)或有限步運(yùn)算，對(duì)無窮過程進(jìn)行截?cái)啵@樣產(chǎn)生的誤差成為截?cái)嗾`差。例如，求一個(gè)級(jí)數(shù)的和或無窮序列的極限時(shí)，取有限項(xiàng)作為它們的近似，從而產(chǎn)生了截?cái)嗾`差。

2.1.3 其他計(jì)算誤差

由計(jì)算格式、迭代格式、運(yùn)算方式和運(yùn)算次數(shù)等的不同引起的誤差。在運(yùn)算過程中，應(yīng)盡量減少迭代次數(shù)，計(jì)算次數(shù)的增加，會(huì)使誤差得到累積和傳播。

2.2 模型誤差

2.2.1 單元離散誤差

由離散單元的幾何形狀產(chǎn)生的誤差，以直邊代替曲邊、以平面代替曲面。使有限元結(jié)果和理論結(jié)果之間存在誤差，當(dāng)離散的單元個(gè)數(shù)越多，計(jì)算精度就越高。但離散的單元個(gè)數(shù)越多，計(jì)算量就越大。

在求解整機(jī)時(shí)，處理接觸問題是計(jì)算的核心，則需要增加接觸域的單元密度，來提高接觸域的計(jì)算精度。必要時(shí)，還需要在迭代過程中，對(duì)接觸域的局部重新自適應(yīng)劃分網(wǎng)格，來提高局部的計(jì)算精度。

2.2.2 形函數(shù)誤差

形函數(shù)誤差是由位移模式的階次低于實(shí)際形函數(shù)多項(xiàng)式的階次而產(chǎn)生的。如果不滿足形函數(shù)的收斂準(zhǔn)則，其誤差不會(huì)因單元數(shù)目的增多和網(wǎng)格的細(xì)化而消失。

N為單元形函數(shù)；δe為已求解出的節(jié)點(diǎn)位移。

則其誤差為：

整個(gè)結(jié)構(gòu)的能量誤差模型為：

3 實(shí)驗(yàn)驗(yàn)證及分析

3.1 實(shí)驗(yàn)裝置

實(shí)驗(yàn)裝置如圖1所示，共用8個(gè)位移傳感器檢測(cè)變形及1個(gè)壓電陶瓷力傳感器檢測(cè)螺栓預(yù)緊力。預(yù)緊力分為10次遞增施加，每次2 000N。然后分10次遞減卸載，再次檢測(cè)其變形。當(dāng)預(yù)緊力較小時(shí)，加載過程平穩(wěn)，比較準(zhǔn)確。但預(yù)緊力較大時(shí)，螺紋間的摩擦力增大，加載過程有振動(dòng)，不容易控制載荷的大小。經(jīng)過多次實(shí)驗(yàn)發(fā)現(xiàn)，2 000N，4 000N和6 000N的試驗(yàn)結(jié)果比較準(zhǔn)確。

圖1 實(shí)驗(yàn)裝置

3.2 理論分析結(jié)果和實(shí)驗(yàn)結(jié)果的比較

在螺栓預(yù)緊力為2 000N和4 000N時(shí)，法向接觸面壓小于2.5MPa，零件本體的變形較小，接觸表面層表現(xiàn)出明顯的非線性特性。因此，對(duì)2 000 N和4 000N時(shí)的分析結(jié)果和實(shí)驗(yàn)結(jié)果進(jìn)行比較。

通過實(shí)驗(yàn)發(fā)現(xiàn)，在預(yù)緊力為2 000N和4 000N時(shí)，小圓柱體上端面上大多數(shù)節(jié)點(diǎn)的位移都為負(fù)值，即小圓柱整體向下位移，與耦合分析的結(jié)果基本一致。

用傳感器檢測(cè)小圓柱體上端面上A，B和C3個(gè)點(diǎn)的位移，然后和理論分析結(jié)果進(jìn)行比較，如圖2所示。

當(dāng)預(yù)緊力從2 000 N增加到4 000N時(shí)，主要比較A，B和C 3個(gè)點(diǎn)的位移增量值，如表1所示。經(jīng)過多次實(shí)驗(yàn)檢測(cè)，取其平均值。

圖2 測(cè)量點(diǎn)的位置

從表1可以看出，耦合分析結(jié)果和實(shí)驗(yàn)測(cè)量結(jié)果非常接近，即零件的變形趨勢(shì)大致相同，但在數(shù)值上有一定的誤差。其原因有以下幾點(diǎn)：

a.由于磨削加工方法和加工精度的不同，零件的表面完整性會(huì)有差異。在相互配對(duì)構(gòu)成結(jié)合面時(shí)，紋理方向的不同也會(huì)影響實(shí)驗(yàn)結(jié)果。因此，結(jié)合面的基礎(chǔ)特性參數(shù)和零件真實(shí)的接觸特性之間存在誤差。沒有考慮接觸域的真實(shí)微觀形狀，按照理論設(shè)計(jì)尺寸建立有限元模型，算法依賴于結(jié)合面的基礎(chǔ)特性參數(shù)，分析結(jié)果是一個(gè)預(yù)測(cè)值。因此，和實(shí)驗(yàn)測(cè)量值之間存在誤差。

b.結(jié)合面的基礎(chǔ)特性參數(shù)基于單位面積，但在加工后零件存在形位誤差，且尺寸越大，誤差就越大。建立有限元模型時(shí)，也沒有考慮零件的形位誤差。因此，理論分析和實(shí)驗(yàn)測(cè)量值之間存在誤差。

c.在實(shí)驗(yàn)時(shí)，由于儀器的測(cè)量精度有限，而且測(cè)量的準(zhǔn)確性也依賴于實(shí)驗(yàn)裝置是否可靠性，以及測(cè)量水平的高低。因此，測(cè)量結(jié)果可能會(huì)不準(zhǔn)確。

表1 耦合分析結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果

4 結(jié)束語

將接觸表面層的接觸面壓與接觸變形非線性關(guān)系，代入到接觸問題的控制方程中，最終讓接觸域的穿透量等于接觸表面層的相對(duì)位移。并討論了實(shí)時(shí)耦合處理和最終耦合處理兩種方法及其適用場(chǎng)合。并且從計(jì)算誤差和模型誤差兩個(gè)方面分析了有限元法的誤差。最后用一個(gè)實(shí)例對(duì)耦合分析軟件進(jìn)行測(cè)試，并用實(shí)驗(yàn)進(jìn)行了多次檢測(cè)，驗(yàn)證了耦合分析結(jié)果與實(shí)驗(yàn)結(jié)果基本一致。因此，提出的機(jī)械系統(tǒng)有限元耦合分析算法是正確的。

［1］ 黃玉美，張廣鵬，高 峰.虛擬樣機(jī)整機(jī)結(jié)構(gòu)特性邊界元仿真［M］.北京：機(jī)械工業(yè)出版社，2004.

［2］ 董獻(xiàn)國.導(dǎo)軌結(jié)合部的變形解析、加工誤差綜合及其結(jié)構(gòu)優(yōu)化設(shè)計(jì)［D］.西安：陜西機(jī)械學(xué)院，1989.

［3］ 王曉春，孔祥安.接觸力學(xué)及其計(jì)算方法［J］.西南交通大學(xué)學(xué)報(bào)，1996，31（3）：230－233.

［4］ 溫衛(wèi)東，高德平.接觸問題數(shù)值分析方法的研究現(xiàn)狀與發(fā)展［J］.南京航空航天大學(xué)學(xué)報(bào)，1994，26（5）：664－667.

［5］ 劉曉平，徐燕申，宋健偉，等.機(jī)械結(jié)構(gòu)結(jié)合面阻尼特征參數(shù)識(shí)別的研究［J］.振動(dòng)與沖擊，1995，14（3）：61－65.

［6］ Ohtakake K，Oden J T，Kikuchi N.Analysis of certain unilateral problem in von karman plate theory by a Penalty method－Part I：a variational primciple with penalty［J］.Computer ＆ Methods in Applied Mechan－ics and Engineering，1980，24（2）：187－213.