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

==概述==

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

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

iSolver介紹視頻：

http://www.yqgqt.org.cn/college/video/c12884

==第20篇：UEL用戶子程序開發步驟==

用戶子程序主要是將用戶特定的材料本構模型和單元算法等公式編寫為計算機語言表示的公式，并實現和商軟求解器之間的交互迭代。

常用的商業有限元軟件都提供了用戶自定義子程序的功能，且一般都是Fortran語言開發，Fortran是上世紀70年代的語言，相對現代化的流行語言編寫，格式要求非常嚴格，編譯調試都比較繁瑣，使得開發效率低下，而且接口限制較多，除了商軟提供的功能外用戶基本沒法改動，靈活性較差。由于用戶子程序很多都涉及復雜的公式編寫，用戶除了需要扎實的理論基礎外，還需要較強的能將公式表達為Fortran語言的編程能力，這對非計算機專業出身的人來說往往在浪費了很多額外精力，使得很多理論高手都對用戶子程序望而卻步，難以入門。

在實際工作中，很多工程師用Matlab來編寫和推導公式，Matlab被認為是市面上最接近草稿紙上推導公式的一款軟件了，而且有限元在數值層面上的計算其實就是矩陣運算，所以Matlab這種數據按矩陣來組織非常適合用來開發有限元相關的程序。而現在市面上還沒有采用Matlab來開發商軟子程序的案例。iSolver是市面上第一款基于Matlab來開發商軟用戶子程序的軟件工具，支持用Matlab編寫和調試用戶子程序。iSolver子程序的接口完全按照Abaqus的標準實現，而Abaqus的子程序接口在近幾年內已經基本不再變化了，同樣的，雖然iSolver在不斷發展，但iSolver子程序接口將維持不變，所有在iSolver上編寫的算法子程序都只要維護自己的算法部分就行，而不是維護整個有限元求解的整個過程。

前面第八、九篇介紹了UMAT用戶自定義材料的開發，這里將介紹UEL用戶自定義單元的開發，本文首先簡單的討論了UEL的一般含義，并詳細的介紹了基于Fortran和Matlab兩種方式的UEL的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。具體的開發過程可以參考下面我們的Step by Step的錄像，包括了整個的有限元基礎理論和我們對Abaqus中的單元理解：

http://www.yqgqt.org.cn/college/video/c14948

深入淺出有限元：基礎理論->Abaqus操作->matlab編程

                                             

1.1 UEL的關鍵輸入輸出參數

UEL網上資料很多，大家可以很容易查看，但大部分資料都只是簡單提供UEL算例，這里我們列出了UEL接口的關鍵輸入輸出參數，如下表所示：

1.2 基于Fortran的Abaqus的UEL的開發步驟

1.2.1 在inp文件中定義UEL

Abaqus中只能通過修改手動inp文件完成用戶UEL的定義，通常包含以下關鍵字及相應屬性，如圖所示：

• *User element關鍵字，用于定義單元基本信息，包括以下屬性：

  nodes為單元節點數；

  type為單元類型，U表示為用戶自定義單元，1001表示單元類型編號，一般從1000開始；

  properties為浮點數屬性個數，一般用于輸入單元對應的section數據和材料數據；

  coordinates為坐標維數；

  variables為狀態量個數，一般求解計算過程中需要在迭代步之間傳遞的變量；

  第二行開始定義節點激活的自由度；

• *Uel property關鍵字，與*User element關鍵字中的properties對應，第二行開始定義具體的數據。

1.2.2 編寫

使用任意編輯器編寫.for文件，推薦使用Visual Studio Code，微軟開源的輕量化代碼編輯器，配置靈活高效：

1.2.3 編譯（可選）

Abaqus沒有自帶Fortran編譯器，所以用戶需要自己去安裝Fortran編譯器和Visual Studio Build Tools，并配置相應環境。具體配置過程與UMAT一致，可以查看我們關于環境配置的視頻：

http://www.yqgqt.org.cn/college/video/c13034?chapter=1

在環境配置完成之后，打開命令提示框，輸入命令Abaqus make Library=XXX.for，即開始編譯，編譯過程中的警告和錯誤都會打印在命令提示框內。

1.2.4 運行

運行有兩種方法，第一種就是在命令提示框中輸入Abaqus job=XXX user=XXX.for，如下圖所示。

第二種就是在Abaqus中創建基于inp文件的任務，然后選擇對應的用戶子程序for文件，在任務管理器中提交運行，如圖所示。

至此，基于Fortran的UEL開發流程已經完成，但結果的正確性還需要更加細致的驗證，為更方便的查找問題，建議先采用單個單元調試UEL，在確保單個單元正確后再將UEL用于實際問題。

1.2.5 調試（可選）

如果想要知道代碼的運行結果是否和預期的一致，一種笨辦法是用print打印到log文件中，高效的方法是采用斷點調試的方法進行運行中的調試。

Abaqus支持命令行調試，不過命令行反復運行也比較繁瑣，用戶也可選擇用一鍵調試Abaqus的用戶子程序的DUS插件工具。DUS（Debug User Subroutine）是集成在ABAQUS/CAE中的一個插件，能夠一鍵啟動用戶配置的用戶子程序開發平臺（如Visual Studio 2008等），并進入對用戶定義子程序的單步調試模式。

有興趣的可到下面網頁下載使用。

Abaqus用戶子程序調試插件：

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

1.3 基于Matlab的iSolver的UEL開發步驟

基于Matlab的Abaqus的UEL具體開發步驟和Abaqus類似，只不過某些步驟需要用到自研有限元求解器開發平臺iSolver。

1.3.1 在inp文件中定義UEL

與Abaqus相應的操作一致，如圖所示：

1.3.2 編寫

在Matlab中創建并編寫U1001.m的文件，放入Abaqus工作目錄下。該文件只包括一個U1001函數，接口和Abaqus的接口參數完全一致，功能也是計算應力應變關系和當前應力狀態等，相對Fortran，利用Matlab可以更容易的編寫計算公式，同時可以利用Matlab在矩陣計算中各種強大功能和算法庫。因為Abaqus的UEL接口和計算功能各個版本相對固定，這個matlab的UEL接口參數也相對固定，不會因為iSolver的版本不同而重新修改接口。

1.3.3 編譯（無）

由于matlab是腳本語言，不需要編譯。

1.3.4 調試（可選）

在Abaqus菜單欄的Plug-ins里選擇iSolver插件的菜單。

點擊iSolver->Engine，按照下圖所示，在功能項Use Solver中選擇iSolver，在Source Type里面選擇Matlab，勾選Debug。點擊Submit進行調試運行。

程序會自動打開matlab并加載U1001.m文件，手動打上斷點

點擊在Debug菜單下的Run U1001運行。

程序將在斷點處停止，且將鼠標移動到需要調試查看的參數上，能夠查看到對應的值。

按F10可以進行單步調試。

1.3.5 運行

在上述步驟的基礎上去掉勾選Debug選項，點擊Submit運行計算，此時將采用iSolver求解器聯合U1001.m進行求解分析，運行完畢點擊Result在Abaqus中查看結果。

1.4 算例

具體的殼理論和在iSolver中實現的UEL算例可以參考下面的視頻：

http://www.yqgqt.org.cn/college/video/c14948

深入淺出有限元：基礎理論->Abaqus操作->matlab編程

 

1.5 總結

本文首先簡單的討論了UEL的一般含義，并詳細的介紹了基于Fortran和Matlab兩種方式的UEL的開發步驟，對比發現開發步驟基本相同，但Matlab更加高效和靈活。同時，由于iSolver基本單元類型和Abaqus算法完全一致，可以發現同一個算例驗證兩者分析結果完全一致，從而證明基于Matlab的UEL的流程和結果的正確性。

UEL的開發一方面要有扎實的公式推導能力，另一方面需要基礎的編程能力和開發工具應用水平，后者不是重點，但往往浪費了大家很多的精力，善用工具方能提高效率，基于Fortran和Matlab兩種方式的UEL的開發步驟和開發工具如下表：

 

如果有任何其它疑問或者項目合作意向，也歡迎聯系我們：

SnowWave02 From www.yqgqt.org.cn

email: snowwave02@qq.com

以往的系列文章：

第一篇：S4殼單元剛度矩陣研究。介紹Abaqus的S4剛度矩陣在普通厚殼理論上的修正。

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

第二篇：S4殼單元質量矩陣研究。介紹Abaqus的S4和Nastran的Quad4單元的質量矩陣。

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

第三篇：S4殼單元的剪切自鎖和沙漏控制。介紹Abaqus的S4單元如何來消除剪切自鎖以及S4R如何來抑制沙漏的。

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

第四篇：非線性問題的求解。介紹Abaqus在非線性分析中采用的數值計算的求解方法。

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

第五篇：單元正確性驗證。介紹有限元單元正確性的驗證方法，通過多個實例比較自研結構求解器程序iSolver與Abaqus的分析結果，從而說明整個正確性驗證的過程和iSolver結果的正確性。

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

第六篇：General梁單元的剛度矩陣。介紹梁單元的基礎理論和Abaqus中General梁單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus梁單元完全一致的剛度矩陣。

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

第七篇：C3D8六面體單元的剛度矩陣。介紹六面體單元的基礎理論和Abaqus中C3D8R六面體單元的剛度矩陣的修正方式，采用這些修正方式可以得到和Abaqus六面體單元完全一致的剛度矩陣。

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

第八篇：UMAT用戶子程序開發步驟。介紹基于Fortran和Matlab兩種方式的Abaqus的UMAT的開發步驟，對比發現開發步驟基本相同，同時采用Matlab更加高效和靈活。

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

第九篇：編寫線性UMAT Step By Step。介紹基于Matlab線性零基礎，從零開始Step by Step的UMAT的編寫和調試方法，幫助初學者UMAT入門。

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

第十篇：耦合約束（Coupling constraints）的研究。介紹Abaqus中耦合約束的原理，并使用兩個簡單算例加以驗證。

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

第十一篇：自主CAE開發實戰經驗第一階段總結。介紹了iSolver開發以來的階段性總結，從整體角度上介紹一下自主CAE的一些實戰經驗，包括開發時間預估、框架設計、編程語言選擇、測試、未來發展方向等。

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

第十二篇：幾何梁單元的剛度矩陣。研究了Abaqus中幾何梁的B31單元的剛度矩陣的求解方式，以L梁為例，介紹General梁用到的面積、慣性矩、扭轉常數等參數在幾何梁中是如何通過幾何形狀求得的，根據這些參數，可以得到和Abaqus完全一致的剛度矩陣，從而對只有幾何梁組成的任意模型一般都能得到Abaqus完全一致的分析結果，并用一個簡單的算例驗證了該想法。

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

第十三篇：顯式和隱式的區別。介紹了顯式和隱式的特點，并給出一個數學算例，分別利用前向歐拉和后向歐拉求解，以求直觀表現顯式和隱式在求解過程中的差異，以及增量步長對求解結果的影響。

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

第十四篇：殼的應力方向。簡單介紹了一下數學上張量和Abaqus中殼的應力方向，并說明Abaqus這么選取的意義，最后通過自編程序iSolver來驗證殼的應力方向的正確性。

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

第十五篇：殼的剪切應力。介紹了殼單元中實際的和板殼近似理論中的剪切應力，也簡單猜測了一下Abaqus的內部實現流程，最后通過一個算例來驗算Abaqus中的真實的剪切應力。

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

第十六篇：Part、Instance與Assembly。介紹了Part、Instance與Assembly三者之間的關系，分析了Instance的網格形成原理，并猜測Abaqus的內部組裝實現流程，隨后針對某手機整機多part算例，通過自編程序iSolver的結果比對驗證我們的猜想。

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

第十七篇：幾何非線性的物理含義。介紹了幾何非線性的簡單的物理含義，并通過幾何非線性的懸臂梁Abaqus和iSolver的小應變情況的結果，從直觀上理解幾何非線性和線性的差異。

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

第十八篇：幾何非線性的應變。首先從位移、變形和應變的區別說起，然后通過一維的簡單例子具體介紹了幾何非線性下的應變的度量方式，并給出了工程應變、 真實應變、Green應變三者一維情況下在數學上的表達方式。

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

第十九篇：Abaqus幾何非線性的設置和后臺。首先介紹了幾何非線性一般的分類，然后詳細說明了Abaqus中幾何非線性的設置方式和常用單元的分類，最后以一個殼單元的簡單算例為對象，可以發現應變理論、Abaqus和iSolver三者在線性、小應變幾何非線性和大應變幾何非線性三種情況下都完全一致，從而驗證Abaqus幾何非線性后臺采用的應變和我們的預想一致。

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