機械系統靜態特性的有限元耦合分析研究

劉 潔

（延安大學化學與化工學院，陜西 延安716000）

0 引言

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

1 有限元耦合分析

1.1 建立接觸系統的求解方程

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

U1為系統Ⅰ非接觸邊界節點的位移列陣；U2為系統Ⅰ接觸邊界節點的位移列陣；F1為系統Ⅰ非接觸邊界節點的力邊界條件列陣；F2為系統Ⅰ接觸邊界節點的接觸力列陣。

U3為系統Ⅱ接觸邊界節點的位移列陣；U4為系統Ⅱ非接觸邊界節點的位移列陣；F3為系統Ⅱ接觸邊界節點的力邊界條件列陣；F4為系統Ⅱ非接觸邊界節點的接觸力列陣。

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

1.2 結合部子結構的罰函數法有限元

將接觸區域的非嵌入條件以及其他附加條件作為懲罰項引入接觸域子結構的泛函中，將接觸問題轉化為關于接觸剛度的優化問題［2－4］。接觸剛度的物理意義是結合面抵抗變形的能力，相當于在接觸點對之間嵌入1個法向、2個切向壓縮彈簧，來抵抗節點間法向和切向的相對位移。

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

建立結合部子結構的泛函：

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

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

當泛函∏C取極值時，即

得到：

整理合并后得到接觸域子結構的控制方程為：

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

罰函數法的優點是算法簡單，易于編程［6］。但接觸條件只是近似被滿足，在接觸域存在穿透，最大相對位移量由給定的常量ε來控制。當選用的罰因子δk過大時，此時的接觸剛度會很大，從而引起計算結果的病態。

1.3 對接觸域的精確處理

如果對計算精度的要求比較高時，需要檢驗罰函數法的計算結果是否合理，并做進一步精確處理。

設第k個接觸點對之間由接觸表面層變形產生的法向位移為λ2k，切向相對位移為λ1k，λ3k，則

合并寫成列陣得到：

則第k個接觸點對之間的接觸面壓為：

將非線性接觸面壓等效為節點載荷得：

［N］為單元形函數矩陣。

將結合部的相對變形｛λk｝和接觸力｛FC｝進行集成裝配，得到總的列陣｛λ｝和｛FC｝。

聯立式（6）和式（12），使用增量迭代法求解，最終讓接觸表面層的最大相對位移量等于接觸表面層的相對變形。當滿足收斂條件后，由式（12）得到接觸節點的接觸力FC列陣，將FC帶入式（1）和式（2）即可求解機械系統Ⅰ與Ⅱ。

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

機械系統Ⅰ和Ⅱ都應滿足整體附加平衡條件：

因結合部相對位移的數量級一般小于等于微米級，而接觸區域的零件宏觀變形和本體結構及外載荷相關，進行耦合時有以下2種迭代方案。

a.實時處理。當接觸區域的零件宏觀變形和結合部的相對變形處于同一個數量級（或者相差較小）時，則在每一次迭代過程中均代入結合面的基礎特性參數，對接觸域進行實時修正。

b.最終處理。當接觸區域的零件宏觀變形和結合部的相對變形相差幾個數量級時，先不代入結合面的基礎特性參數，即先用罰函數法求出大致的接觸域后，再代入結合面的基礎特性參數做準確處理，對接觸域進行精確的修正。

2 誤差分析

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

2.1 計算誤差

2.1.1 舍入誤差

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

2.1.2 截斷誤差

由于實際運算只能完成有限項或有限步運算，對無窮過程進行截斷，這樣產生的誤差成為截斷誤差。例如，求一個級數的和或無窮序列的極限時，取有限項作為它們的近似，從而產生了截斷誤差。

2.1.3 其他計算誤差

由計算格式、迭代格式、運算方式和運算次數等的不同引起的誤差。在運算過程中，應盡量減少迭代次數，計算次數的增加，會使誤差得到累積和傳播。

2.2 模型誤差

2.2.1 單元離散誤差

由離散單元的幾何形狀產生的誤差，以直邊代替曲邊、以平面代替曲面。使有限元結果和理論結果之間存在誤差，當離散的單元個數越多，計算精度就越高。但離散的單元個數越多，計算量就越大。

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

2.2.2 形函數誤差

形函數誤差是由位移模式的階次低于實際形函數多項式的階次而產生的。如果不滿足形函數的收斂準則，其誤差不會因單元數目的增多和網格的細化而消失。

N為單元形函數；δe為已求解出的節點位移。

則其誤差為：

整個結構的能量誤差模型為：

3 實驗驗證及分析

3.1 實驗裝置

實驗裝置如圖1所示，共用8個位移傳感器檢測變形及1個壓電陶瓷力傳感器檢測螺栓預緊力。預緊力分為10次遞增施加，每次2 000N。然后分10次遞減卸載，再次檢測其變形。當預緊力較小時，加載過程平穩，比較準確。但預緊力較大時，螺紋間的摩擦力增大，加載過程有振動，不容易控制載荷的大小。經過多次實驗發現，2 000N，4 000N和6 000N的試驗結果比較準確。

圖1 實驗裝置

3.2 理論分析結果和實驗結果的比較

在螺栓預緊力為2 000N和4 000N時，法向接觸面壓小于2.5MPa，零件本體的變形較小，接觸表面層表現出明顯的非線性特性。因此，對2 000 N和4 000N時的分析結果和實驗結果進行比較。

通過實驗發現，在預緊力為2 000N和4 000N時，小圓柱體上端面上大多數節點的位移都為負值，即小圓柱整體向下位移，與耦合分析的結果基本一致。

用傳感器檢測小圓柱體上端面上A，B和C3個點的位移，然后和理論分析結果進行比較，如圖2所示。

當預緊力從2 000 N增加到4 000N時，主要比較A，B和C 3個點的位移增量值，如表1所示。經過多次實驗檢測，取其平均值。

圖2 測量點的位置

從表1可以看出，耦合分析結果和實驗測量結果非常接近，即零件的變形趨勢大致相同，但在數值上有一定的誤差。其原因有以下幾點：

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

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

c.在實驗時，由于儀器的測量精度有限，而且測量的準確性也依賴于實驗裝置是否可靠性，以及測量水平的高低。因此，測量結果可能會不準確。

表1 耦合分析結果與實驗測量結果

4 結束語

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

［1］ 黃玉美，張廣鵬，高 峰.虛擬樣機整機結構特性邊界元仿真［M］.北京：機械工業出版社，2004.

［2］ 董獻國.導軌結合部的變形解析、加工誤差綜合及其結構優化設計［D］.西安：陜西機械學院，1989.

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

［4］ 溫衛東，高德平.接觸問題數值分析方法的研究現狀與發展［J］.南京航空航天大學學報，1994，26（5）：664－667.

［5］ 劉曉平，徐燕申，宋健偉，等.機械結構結合面阻尼特征參數識別的研究［J］.振動與沖擊，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.