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

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了。

我們將用統一的公式來求解這兩類關系，同時也從軟件實現層面說明一下針對這兩類情況的各自差異。分幾篇文章來介紹約束關系，上一章我們介紹了基于點或者面之間的約束關系的形式是統一的關系，這章我們就將以此具體推導這種統一的約束關系在有限元中如何采用Lagrange因子法求解。

3 Lagrange因子法求解約束方程

沒加約束前的虛功原理表示為：

4 各個軟件的約束關系的具體實現方式

約束關系雖然原理都可以用Lagrange因子法來求，但軟件的實現往往更復雜，實際編程有具體的兩種實現方式：

（1） 既然上一章說了約束關系可以類比普通的有限單元，那么可以把約束關系按單元來實現。

（2） 約束關系也就是Master節點自由度和Slave節點自由度之間的系數，直接把這個約束關系寫下來放到全局剛度陣就行。

4.1 Nastran的實現方式

Nastran理論手冊中說明了一開始版本兩者都有，但后面更多推薦用單元方式來實現，而且默認情況下Nastran的MPC等約束關系相關關鍵詞后面第一個參數EID是Element Number，如下圖，所以猜測采用的單元形式。

4.2 Abaqus的實現方式

Abaqus的關鍵詞中inp的MPC中已經消除了這個單元號，如下圖，Abaqus如果是按單元來實現的話應該會在inp關鍵詞中反應，譬如質量點，雖然Abaqus/CAE上是inertia設置，但inp依然是Element關鍵詞，也會設置單元類型為MASS，但MPC猜測沒有這么做，而是把Master點和Slave點的關系寫到了關系矩陣中，在組裝全局剛度陣時將關系式單獨寫到全局剛度陣中。

4.3 iSolver的實現方式

iSolver實現時一開始是按第二種方法，也就是Abaqus的方式來做的，但發現這種方法當遇到復雜問題時改動比較大，無論哪種方法，Master節點自由度和Slave之間的關系都是局部節點的系數，這個系數在全局節點下都需要進一步轉換，而第二種方法就需要再自己轉換一次，當遇到復雜問題譬如兩個約束關系嵌套（一個約束關系的Master是另一個約束關系的Slave）、或者一個Slave被兩個Master用到時這個轉換的編程復雜度也會上去，后面我們對于部分約束關系改為了單元形式，因為此時就只需要在整體組裝時一次性的將局部節點自由度關系轉換到全局關系了。也許可能是我們編程水平不夠才沒有做到和Abaqus一樣的方法，具體Abaqus怎么做的要是哪位大神知道也歡迎聯系我們。

如果你自己編程序，不管采用哪種實現方式，我們覺得除了要實現這個功能外，還要考慮你的程序的兼容性，不破壞你原有的架構。

軟件具體實現的問題就不再詳述了，我們再接著討論兩個Lagrange因子求解約束方程兩個基本的容易搞混的問題：

5 約束方程的個數

約束關系具體是各個節點自由度之間的關系，但為了簡單起見，我們只認為是節點之間的約束關系，拋開自由度個數，也就是無論節點是3個還是6個，只要還是這個節點的約束方程，個數只認為是1個。

那么每種約束關系（每個節點下的具體的每個自由度約束對應的一個Lagrange因子）有多少個約束方程？

（1） 點之間的約束關系

RBE2每個Slave節點的位移都必須和Master節點相等，也就是每個Slave節點對應1個約束關系，所以RBE2的約束關系數目=Slave節點數，對每一個Slave節點自由度，都應該有一個Lagrange因子，該Lagrange因子可以認為綁定在Slave節點上的物理量。而RBE3中因為所有的Slave和一個Master節點組成了一個約束關系，所以只有一個真正的約束關系，所以RBE3的約束關系數目=Master節點數（一般是1個），Lagrange因子可以認為綁定在Master節點。

（2）面之間的約束關系

由面之間的約束關系方程，可知對每個Slave節點，都有一個約束方程，因此， Contact或者Tie等面之間的約束方程個數=Slave節點數目，和RBE2一致。Lagrange因子可以認為綁定在Slave節點。

6 Lagrange因子的物理含義

由上述方程從形式上來看，可知，Lagrange因子類比普通單元中的應力， 而h類比普通單元中的應變，但還是很多人無法理解Lagrange因子到底是什么，我們理解也是實際的物理量。分兩種情況

6.1 點之間約束的Lagrange因子含義

對于點之間的約束關系的h為Master點和Slave點的位移差（簡單起見，我們不考慮轉動），單位是長度

6.2 面之間約束的Lagrange因子含義

對于面之間的約束關系的h為Master面和Slave面的方向距離，單位是長度

約束關系的能量項為：

6.3 點之間約束的Lagrange因子含義證明

為證明點約束的Lagrange因子是集中力，我們取一個體單元43.75X10X5及一個RP點，RP點添加x方向集中力載荷1e9。

RP點和體單元采用KCoupling綁定。

6.3.1 單節點綁定

第一次僅綁定8號節點。

采用iSolver求解Lagrange因子得到如下（此時沒有打開幾何非線性）：

可以發現8號節點Lagrange因子為-1e9、1.95e-3、1.95e-3，相對e9來說，yz方向的Lagrange因子可以略為0，得到該節點正好和右邊的x方向的力平衡，和前面面之間約束關系類似，KCoupling的Lagrange因子為Slave節點處Slave節點對Master節點的作用力，在圖上表示如下。方向或者正負號完全由約束方程的正負決定，iSolver是這種負號是因為取的約束方程是h=Us-Um，而不是h=Um-Us。

而且這個作用力在當前時刻猜測不是沿著Slave節點指向Master節點，因為此時iSolver計算出來的節點8的位移在yz方向要明顯小于RP點位移，所以，8節點指向RP點已經不是沿x方向了，也就是Slave和Master節點之間不是類似桿只承受拉壓，而是類似梁可以承受切向力。

6.3.2 多節點綁定

Slave采用7號和8號節點綁定。如下圖所示：

采用iSolver求解Lagrange因子得到：

可以發現7號節點的z正好和8號節點z方向作用力抵消，而8號節點的x方向的作用力正好和外力抵消。

其實KCoupling綁定的所有Slave節點Lagrange因子的和肯定和外力相同，這點可以直接用約束關系推導求出，有興趣的可以試一下。

6.3.3 問題

上述點之間約束關系的Lagrange是Slave節點和Master節點之間的作用力怎么用Abaqus來驗證一直沒有找到類似接觸的CPRESS合適的對應物理量，要是哪位大神知道也請告知一下。

6.4 面之間約束的Lagrange因子含義證明

    為證明面約束的Lagrange因子是單位面積上的接觸力，我們取兩個5X1X0.1的長方體，設置兩個面緊貼在一起，設置不滑動，下方Master體單元的四個點固支，上方Slave單元的四個點加法向1e10的壓力。

6.4.1 理論值

顯然，理論上的面壓力=1e10*4/(5X1)=8e9。

6.4.2 iSolver解

在iSolver中，我們求取Lagrange因子采用和Slave節點相同的自由度所在全局中的位置。查詢單元，可知Slave單元為1號單元，由5,6,8,7,1,2,4,3組成，而Slave節點為8,7,4,3。每個節點三個自由度，即Lagrange因子所在自由度為7-12，和19-24。

采用iSolver求解器求解Lagrange因子，打印出Lagrange的值：

可發現第9和12的Lagrange因子的值為-8e9，和理論一致。由于是負值，所以指向Z-方向，表示Slave節點處Slave節點對Master面的壓力。

6.4.3 Abaqus解

將iSolver得到的模型輸出到Abaqus進行求解，暫時沒找到Abaqus直接的Lagrange因子打印方法，但既然我們認為Lagrange就是面力，可以查看Abaqus接觸面力輸出。

Abaqus接觸面壓力為CPRESS，查看CPRESS如下圖所示，為8e9，猜測CPRESS是負的Lagrange因子，表示Slave節點處Master對Slave節點的面壓力。

同時，我們也可以查看Abaqus的單位面上的切向摩擦力，即CSHEAR1/2，如下所示：

顯然，和iSolver的Lagrange因子正好差個負號，猜測也是取了負的Lagrange因子。

7 以往的系列文章

7.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

7.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

7.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

7.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

7.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

 第四十五篇：約束關系（1）-統一形式

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