（原創，轉載請注明出處）

1 概述

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1） 基礎理論

（2） 商軟操作

（3） 自編程序

三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。

有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。

一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元軟件iSolver，通過自研CAE軟件和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。iSolver包括完整的前后處理和有限元求解器，功能如下，有興趣可直接在下面網址下載：

百度網盤鏈接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取碼: 6fdf

2 約束關系

我們在系列文章第35章介紹了接觸分析及常用的基本算法。現在Abaqus、LS-DYNA、Ansys等結構商軟都說可以處理復雜的上萬零部件接觸的整車、整機等模型仿真，沒做過實際的這種仿真分析，很好奇，接觸分析算法往往涉及大變形、邊界不連續，只要輸入條件或者算法稍微變化一些，兩個零部件算出來的接觸結果就可能差異很大，更不用說上萬個零部件的接觸結果了，對這種大規模組裝模型的仿真結果不知如何來判斷它的可靠性，像普通的只校核一下材料的應力還是看一下動畫是否和試驗一致？畢竟仿真只有簡單的標準來判斷結果的正確性才能在企業中起到真正輔助設計的作用，如果你恰好做過，不妨也簡單介紹一下你的經驗。對自研CAE軟件開發者來說，自研結構CAE軟件是否真的要和商軟去比拼接觸等復雜算法還是花更多時間在精雕細琢那些常用功能上，這也是開發者需要慎重考慮的問題，而且很多自主CAE軟件連常規線性問題都算的不對，或者都沒法用鼠標穩定的走完那些材料、屬性、邊界、加載等流程，用戶又怎么會相信你能算對接觸這種復雜問題的？

不管怎么樣，從有限元實現的角度來講，如果想做真正實際工程中的接觸分析，那么首先需要去研究約束關系，接觸分析在有限元中也僅是約束關系的一種。有限元中的約束很多場景大家用的是邊界中的簡支、固支等約束，但從更廣泛的角度上講，只要表示一個節點的某個自由度依賴于其它的節點自由度或者取某個特定值，就可以稱為約束關系。只不過對固支、簡支等直接自由度=0，在有限元中直接減縮剛度陣就行，很容易求，但對節點自由度相互依賴的約束關系就比較復雜了。約束關系主要有兩類。

（1） 一類是MPC點之間的約束。Nastran的MPC的靈活度要遠遠超過Abaqus，Nastran的主節點可以選擇123自由度，也可以對每個從節點設置不同的自由度，還能主節點和從節點互相包含，Abaqus更多的是只負責80%的常用應用場景，復雜功能讓你編子程序，但事實上一線仿真工程師又有多少人愿意編子程序呢？這種做法導致雖然Abaqus無論從用戶體驗、非線性還是商業化都比Nastran好很多，但很多線性的工程復雜問題還是沒法替代Nastran。

（2） 另一類是Contact、Tie等的面之間的約束關系。在這方面Abaqus要明顯強于Nastran了。

我們將用統一的公式來求解這兩類關系，同時也從軟件實現層面說明一下針對這兩類情況的各自差異。分幾篇文章來介紹約束關系，本篇是約束關系（1）-統一形式，既然接觸僅是約束關系的一種，那么MPC、Tie、接觸等的求解過程也是很類似的，這里將介紹一下這些約束關系如何表達為統一形式。

3 統一形式的約束關系

在沒有約束關系時，如下圖情況，物體在體外力和面外力作用下變化。

有限元方程按照虛功原理求解，在物理上可解釋能量守恒原理，即在某一個時刻點，假定在外力作用下有個虛擬的位移，那么外力在虛擬位移下做的虛功=內部應變能的變化相同。

虛功原理中的每項都表示各自區域在虛位移下的能量變換

（1） fv是每單位體積內的力，外力fv和位移相乘表示單位體積內的虛功，所以對體積積分

（2） fs是每單位面積上的力，外力fv和位移相乘表示單位面積上的虛功，所以對面積積分，推論就是最后一項應力和應變是單位體積內的內能，所以對體積積分。

當存在約束關系時，在能量中加入約束關系相關的一項：

顯然的單位也是該約束關系所在區域在虛位移下的能量單位。

約束關系在有限元中可以分為兩大類：

3.1 點之間的約束關系

（1）點之間約束關系，最常見的是節點之間的剛性連接，Nastran中稱為RBE2，在Abaqus或者iSolver中稱為Kinematic Coupling，此時可以認為Master節點和Slave節點之間焊死在了一個剛性無窮大的直桿上。在實際情況中，Slave節點之間沒有相對位移，但由于計算時很多時候默認為小位移，反而導致Slave節點之間是有相對位移的。

此時：

其中，

為Master節點的旋轉向量，為從Slave節點到Master節點的向量。

這種節點之間的約束關系還有一種常見情況就是節點間的分布耦合連接RBE3，將施加在Master節點上的力和力矩按照加權系數分配到Slave節點上，從而實現載荷在單元間的傳遞，此時就避免了RBE2太剛的問題，典型應用場景譬如模擬圓管對壁面的壓力作用。

分布耦合連接是將Master節點的力和力矩按某種規則分配到Slave節點，保證力和力矩分別相等，即：

上述對六自由度的Master節點來說，一共只有六個方程，但每個Slave節點都有6個位置量，所以對多個Slave節點情況解不唯一，Abaqus、iSolver、Nastran等都有各自不同的分配原則，一般都是假定Slave節點不再存在Master節點分配過來的彎矩Ms，同時，Fs等比例分配。

對這種點的約束關系，虛功原理中增加的能量項表示為：

注意，因為是點之間的約束，所以不需要對面或者體積分，類似質量點對質量陣的貢獻或者集中力對載荷向量的貢獻時也不需要積分一樣。

3.2 面之間的約束關系

面之間的約束關系，最常見的就是Master面和Slave面/點的接觸連接關系。

在有限元軟件中設置完接觸關系對后，對Slave面上的任意一點xs，在Master面上尋找和xs距離最短的點，譬如xm點，此時xs和xm之間的距離為h，顯然，Master面在xm點上的法向n指向xs點（要不然就不是最近鄰點了）。也就是說slave面上任意一點在Master面上都有一個點來對應，這個點我們稱為錨點，錨點的法向指向Slave點是一個必要條件（不是充分條件，因為有可能Master有兩個錨點法向指向同一個xs點，譬如Master是個拐角的情況）。如下圖所示。

此時距離(真實的距離肯定是正的，但我們這邊為了方便取有向距離)

設接觸面法向壓力為p，對硬接觸，當兩個面之間的距離h是0時，壓力p>0。有距離h<0時，沒有接觸關系，此時壓力p為0。

即約束方程為

這是對面上的任意一點。計算機是無法處理面上所有點的問題的，所以也只能劃分成有限單元，得到有限節點之間的關系。

劃分完網格后，上述接觸面上必然也是節點組成，如下：

對每個Slave節點xs，在Master接觸面上尋找對應的錨點，因為Master面由不連續的節點組成，所以這個錨點很多情況都不在Master節點上。譬如上面圖示，此時就是該點所在單元對錨點的插值，同樣滿足這個錨點的法向n(r,s)（也由所在單元的節點在r，s點插值得到）指向xs。

同樣，此時的約束方程為：

這種面面之間的一種特殊情況就是Master面和Slave面/點的粘貼連接關系，表示Slave面/點用膠水粘在Master面上，所以在Nastran中稱為Glue，在Abaqus中稱為Tie，一般用于兩個不同網格之間邊界耦合在一起，譬如下方的橋身和橋墩，實際上Slave面/點是焊死在Master面的，Slave面/點會隨著Master面一起拉伸壓縮移動，但Slave面/點在運動過程中不會脫落，最終反應到剛度陣依然是固定的Slave節點自由度和Master節點自由度之間的約束關系。

對這種面與面/點的約束關系，增加的能量項表示為：

Sc為接觸面。

上式和節點的約束關系相比，可知pms就是Lagrange因子。

在實際三維接觸分析中，接觸力除了法向n的壓力還有兩個切向v1和v2的摩擦力。得到實際接觸力由三個正交的力分量組成：

簡單起見，我們假定不滑動，此時由于接觸關系增加的能量項為：

4 統一形式的約束關系

上述不同的僅是約束關系帶來的積分范圍和約束方程的個數和物理量。去掉所有的積分形式，我們可以統一寫成如下形式:

其中上面的正負不影響結果，僅僅和h的正負有關。

5 更大統一形式的有限元的約束關系

MPC、接觸等都是某部分節點和另一部分節點之間的關系，那么有限元中梁、殼、體普通的單元是什么呢？從更大的統一形式來說，這些普通單元也僅僅是節點之間的約束關系而已，只不過這個關系可能更加復雜點。

以兩節點梁單元為例。

此時，每個節點6個自由度，那么單元剛度陣K為12X12的矩陣。K可以按6x6矩陣分塊，每塊小矩陣都比較類似，那么可以分成4塊：

（1）左上角為第一個節點自由度相關的6X6的小矩陣做研究對象。iSolver軟件中內部打印如下

K的任意一項Kij都是i和j兩個自由度的乘積，其中對角元素為節點1和節點1自身的自由度的乘積，譬如K11就是u和u的乘積，而非對角元素都是耦合項，譬如K15就是u和θy的乘積。

（2）右下角類似，只不過是節點2和節點2自身的關系。

（3）左下角和右上角類似，是節點1和節點2的關系。

所以，可以總結有限元的所有單元，包括MPC、接觸等最終計算的都是這個單元所涉及點兩兩之間的關系，這個關系的最終形式都體現在剛度陣上。有限元簡單來講就是求各個點之間關系的公式，只不過有些公式簡單，譬如兩個點之間用線性彈簧，有些復雜，譬如幾何大變形時兩個節點之間或者接觸時兩個面之間的關系。

6 以往的系列文章

6.1 ========第一階段========

第一篇：S4殼單元剛度矩陣研究。

http://www.yqgqt.org.cn/content/post/338859

第二篇：S4殼單元質量矩陣研究。

http://www.yqgqt.org.cn/content/post/343905

第三篇：S4殼單元的剪切自鎖和沙漏控制。

http://www.yqgqt.org.cn/content/post/350865

第四篇：非線性問題的求解。

http://www.yqgqt.org.cn/content/post/360565

第五篇：單元正確性驗證。

http://www.yqgqt.org.cn/content/post/373743

第六篇：General梁單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/403932

第七篇：C3D8六面體單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/430177

第八篇：UMAT用戶子程序開發步驟。

http://www.yqgqt.org.cn/content/post/432848

第九篇：編寫線性UMAT Step By Step。

http://www.yqgqt.org.cn/content/post/440874

第十篇：耦合約束（Coupling constraints）的研究。

http://www.yqgqt.org.cn/content/post/531029

6.2 ========第二階段========

第十一篇：自主CAE開發實戰經驗第一階段總結。

http://www.yqgqt.org.cn/content/post/532475

第十二篇：幾何梁單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/534362

第十三篇：顯式和隱式的區別。

http://www.yqgqt.org.cn/content/post/537154

第十四篇：殼的應力方向。

http://www.yqgqt.org.cn/content/post/1189260

第十五篇：殼的剪切應力。

http://www.yqgqt.org.cn/content/post/1191641

第十六篇：Part、Instance與Assembly。

http://www.yqgqt.org.cn/content/post/1195061

第十七篇：幾何非線性的物理含義。

http://www.yqgqt.org.cn/content/post/1198459

第十八篇：幾何非線性的應變。

http://www.yqgqt.org.cn/content/post/1201375

第十九篇：Abaqus幾何非線性的設置和后臺。

http://www.yqgqt.org.cn/content/post/1203064

第二十篇：UEL用戶子程序開發步驟。

http://www.yqgqt.org.cn/content/post/1204261

6.3 ========第三階段========

第二十一篇：自主CAE開發實戰經驗第二階段總結。

http://www.yqgqt.org.cn/content/post/1204970

第二十二篇：幾何非線性的剛度矩陣求解。

http://www.yqgqt.org.cn/content/post/1254435

第二十三篇：編寫簡單面內拉伸問題UEL Step By Step

http://www.yqgqt.org.cn/content/post/1256835

第二十四篇：顯式求解Step By Step。

http://www.yqgqt.org.cn/content/post/1261165

第二十五篇：顯式分析的穩定時間增量。

http://www.yqgqt.org.cn/content/post/1263601

第二十六篇：編寫線性VUMAT Step By Step。

http://www.yqgqt.org.cn/content/post/1266640

第二十七篇：Abaqus內部計算和顯示的應變。

http://www.yqgqt.org.cn/content/post/1273788

第二十八篇：幾何非線性的T.L.和U.L.描述方法

http://www.yqgqt.org.cn/content/post/1282956

第二十九篇：幾何非線性的T.L.和U.L.轉換關系

http://www.yqgqt.org.cn/content/post/1286065

第三十篇：諧響應分析原理

http://www.yqgqt.org.cn/content/post/1290151

6.4 ========第四階段========

第三十一篇：自主CAE開發實戰經驗第三階段總結

http://www.yqgqt.org.cn/content/post/1294824

第三十二篇：諧響應分析算法

http://www.yqgqt.org.cn/content/post/1299983

第三十三篇：線性瞬態動力學

http://www.yqgqt.org.cn/content/post/1302074

第三十四篇：非線性瞬態分析

http://www.yqgqt.org.cn/content/post/1787283

第三十五篇：接觸求解算法

http://www.yqgqt.org.cn/content/post/1792869

 第三十六篇：DLOAD用戶子程序開發步驟

http://www.yqgqt.org.cn/content/post/1826803

 第三十七篇：梁單元差異（1）-理論基礎

https://jishulink.com/content/post/1872208

 第三十八篇：梁單元差異（2）-梁截面方向

http://www.yqgqt.org.cn/content/post/1874628

 第三十九篇：梁單元差異（3）-剪力和彎矩

http://www.yqgqt.org.cn/content/post/1876013

 第四十篇：梁單元差異（4）-形心、剪心和偏置

http://www.yqgqt.org.cn/post/1888000

6.5 ========第五階段========

 第四十一篇：自主CAE開發實戰經驗第四階段總結

http://www.yqgqt.org.cn/post/1905963

 第四十二篇：聲學分析(1)-有限元

http://www.yqgqt.org.cn/post/1912491

 第四十三篇：聲學分析(2)-邊界元

http://www.yqgqt.org.cn/post/1923936

 第四十四篇：聲學分析(3)-濕模態

http://www.yqgqt.org.cn/post/1928692