摘 要：ADINA（automatic dynamic incremental nonlinear analysis， ADINA）84作為大型通用有限元程序，具有非常強(qiáng)大的非線性計(jì)算能力，廣泛運(yùn)用于各種工程數(shù)值模擬分析領(lǐng)域，但對其源代碼的解析以及如何運(yùn)用源代碼并對源代碼進(jìn)行二次開發(fā)的研究較少。通過對源代碼深入解析，基于廣義胡克定律編寫材料子程序并編譯了相對應(yīng)的求解器，對比驗(yàn)證了求解器和材料子程序的計(jì)算準(zhǔn)確性，為ADINA源代碼的使用提供了詳細(xì)說明，同時為本構(gòu)模型二次開發(fā)提供參考，為國產(chǎn)專用有限元軟件開發(fā)提供新的思路。

關(guān)鍵詞：ADINA源程序；求解器；幾何非線性；二次開發(fā)

中圖分類號：TB115.1 DOI：10.16375/j.cnki.cn45-1395/t.2024.04.002

0 引言

ADINA（automatic dynamic incremental nonlinear analysis， ADINA）是一款動力非線性有限元分析軟件，且為開源程序[1]，早期國內(nèi)學(xué)者對其進(jìn)行了一定的研究[2-5]。通過幾十年發(fā)展，該軟件已實(shí)現(xiàn)商業(yè)化并運(yùn)用在航空航天、國防軍工、基礎(chǔ)設(shè)施建設(shè)等多個領(lǐng)域。郭富強(qiáng)等[6]基于ADINA對河套灌區(qū)預(yù)制混凝土渠道保溫防凍脹效果進(jìn)行了數(shù)值模擬。貢悅[7]通過ADINA對實(shí)橋試驗(yàn)的拉索減振被動式阻尼器進(jìn)行了數(shù)值模擬。謝建等[8]利用ADINA建立最易失效狀態(tài)下液壓缸模型對導(dǎo)彈的起豎液壓缸進(jìn)行有限元分析。陳志濤[9]針對卵形消化池通過ADINA數(shù)值模擬對該結(jié)構(gòu)系統(tǒng)地震作用下的流固耦合動力特性進(jìn)行了有限元分析，ADINA提供了豐富的用戶接口，使用戶能夠進(jìn)行特殊材料及工況下的分析求解。熊玉春等[10]通過ADINA軟件對3個常規(guī)三軸壓縮模型進(jìn)行材料二次開發(fā)，驗(yàn)證了該軟件具有強(qiáng)大的非線性求解能力；同時工業(yè)軟件的二次開發(fā)同樣重要。方瑞等[11-13]基于UG以及Revit等軟件的二次開發(fā)進(jìn)行了避障仿真研究與安全管理平臺研究等內(nèi)容，在這方面ADINA能力同樣突出。張衛(wèi)東等[14]通過ADINA實(shí)現(xiàn)了黏土廣義本構(gòu)模型二次開發(fā)，為后續(xù)其他彈塑性模型提供參考。汪仁和等[15]改進(jìn)的西原模型生成動態(tài)鏈接庫用于ADINA深井凍結(jié)壁的有限元計(jì)算。由上述文獻(xiàn)可知ADINA具有強(qiáng)大的線性、非線性及多物理場求解功能，同時具有非常強(qiáng)大的二次開發(fā)能力。

但上述文獻(xiàn)都是利用基于ADINA提供的MakeFile文件自動對各個文件鏈接編譯使用，或直接利用ADINA商業(yè)軟件以及ADINA提供的動態(tài)鏈接庫進(jìn)行求解及二次開發(fā)，并且前者編譯環(huán)境老舊，為已經(jīng)停止更新的Compaq Visual Fortran（CVF）開發(fā)環(huán)境。后者與其他商業(yè)有限元軟件一樣不提供開源代碼，很少有直接說明源代碼使用方法及二次開發(fā)的案例。

在自主化趨勢越發(fā)明顯的今天，涌現(xiàn)出多種國產(chǎn)CAE軟件，如NWorks、Simright Simulator、FEPG等，但都比較針對某一特定方向或功能不夠完善。為此，本文解析了ADINA84源代碼實(shí)現(xiàn)流程及二次開發(fā)途徑，基于Window10系統(tǒng)編譯求解器，根據(jù)文獻(xiàn)[16]模型建立輸入文件進(jìn)行大變形分析，以及根據(jù)廣義胡克定律編寫了線彈性材料子程序嵌入求解器，進(jìn)而對8節(jié)點(diǎn)單元模型進(jìn)行求解，驗(yàn)證了所編譯求解器的準(zhǔn)確性以及子程序的可行性，且該源代碼為開源代碼，可以為國產(chǎn)有限元軟件的發(fā)展提供參考。

1 ADINA程序編譯及計(jì)算流程

1.1 ADINA求解器編譯

ADINA源代碼為FORTRAN77語言編寫，與最新的FORTRAN語言版本相比較為老舊，若在如今的電腦系統(tǒng)下成功編譯求解器需要提前安裝相對應(yīng)的FORTRAN編譯器。現(xiàn)在主流的編譯器有以下幾種：Intel Visual Fortran、PGI Fortran、NAG Fortran以及GFortran。其中Intel Visual Fortran優(yōu)化最佳，并且其配套集成開發(fā)環(huán)境最為完善，適合商業(yè)軟件開發(fā)，但容易出現(xiàn)Bug。為避免在Windows環(huán)境下出現(xiàn)不必要的Bug，以及為適應(yīng)研究所需使程序移植便利，采用GFortran開發(fā)環(huán)境，此開發(fā)環(huán)境完全開源的同時不需要額外安裝配套的開發(fā)環(huán)境并且優(yōu)化能力較強(qiáng)，適合學(xué)習(xí)、研究及常規(guī)開發(fā)。

ADINA84源代碼具體編譯過程如下：

Step 1 安裝Minimalist GNU for Windows工具集，其中包含gcc，C++以及Fortran編譯器gfortran。

Step 2 更改賬戶系統(tǒng)環(huán)境變量，添加工具集中g(shù)fortran文件路徑到環(huán)境變量Path路徑中。通過系統(tǒng)命令提示符輸入gfortran-version，若提示Built by MinGW-W64 project，則表明編譯環(huán)境搭建成功。

Step 3 通過編譯器鏈接源代碼所在文件夾的44個源代碼文件，可得到一個自命的.exe的可執(zhí)行文件，整個流程如圖1所示。

1.2 輸入文件

編譯成功的求解器可讀取.txt以及.dat等各種類型的文本文件，但文本中數(shù)據(jù)格式需要滿足ADINA84源代碼的讀取要求；因ADINA84源代碼為FORTRAN77編寫，所以輸入文件格式也是按照“卡片”的形式來讀取信息，其中每一張“卡片”或每一個“卡片組”都是靠卡片名以及卡片所用格式、變量名稱、變量的含義以及注釋來進(jìn)行識別。

在語句格式中，A代表字符變量，I代表整型變量，F(xiàn)或E表示實(shí)型變量。整型變量名以字母I、J、K、L、M、N開頭，字符型以及實(shí)型變量以其他字母開頭。其中整型變量可以由一個“+”或“-”號和數(shù)字組成，由于空格會被讀為0，所以整數(shù)在所留位置中需要從右往左輸入，同樣地，實(shí)型變量在沒有小數(shù)點(diǎn)的情況下也需要靠右輸入。

數(shù)據(jù)結(jié)構(gòu)與現(xiàn)在主流的CAE軟件輸入文件類型類似，例如ABAQUS的.inp文件等，主要包括結(jié)構(gòu)控制卡、載荷控制卡、時間函數(shù)控制卡、節(jié)點(diǎn)控制卡、單元控制卡等。主要的輸入方式（按列號由右往左輸入）、變量名及其意義如表1—表5所示。

1.3 求解器求解總流程

ADINA84源代碼與ADINA81源代碼除部分變量以及二次開發(fā)接口外有許多相似之處[17]，其中計(jì)算主要分為5部分：1）文件讀入、儲存和初始化單元節(jié)點(diǎn)信息、求解狀態(tài)；2）剛度矩陣計(jì)算；3）平衡迭代方式；4）計(jì)算應(yīng)力；5）輸出結(jié)果?？傆?jì)算流程框架見本文二維碼。

2 源代碼二次開發(fā)

2.1 源代碼二次開發(fā)要點(diǎn)

ADINA84源代碼中提供2D和3D實(shí)體單元的子程序接口，其中子程序SUBROUTINE CUSER2針對2D實(shí)體單元，子程序SUBROUTINE CUSER3針對3D實(shí)體單元，且即使自定義的材料不含和溫度有關(guān)的參數(shù)，在實(shí)際建模中也需要對給定材料的溫度和模型的溫度荷載賦值。與ADINA81源代碼不同的是無需定義歷史工作數(shù)組長度，ADINA84二次開發(fā)子程序包括4個程序段，分別對應(yīng)于程序中的KEY1—KEY4。在KEY1程序段，初始化積分點(diǎn)上存儲的變量；在KEY2程序段，對每個STEP的每個DEPS調(diào)用1次，對積分點(diǎn)進(jìn)行應(yīng)力積分并完成歷史變量的更新；在KEY3程序段，計(jì)算彈塑性矩陣；在KEY4程序段，輸出所需的應(yīng)力-應(yīng)變以及位移等結(jié)果，其中自定義子程序代碼需寫在CONTINUE與RETURN之間，公用變量（彈性模量等）可寫在GO TO之前成為全局變量，每個GO TO語句對應(yīng)的KEY程序段之間相互獨(dú)立對應(yīng)不同的計(jì)算子程序，不能隨意刪除，否則容易出現(xiàn)重復(fù)疊加計(jì)算等其他問題。其中材料子程序接口常用變量及解釋見表6。

2.2 子程序編譯思路

1）KEY1初始化ARRAY以及IARRAT兩個數(shù)組。

2）KEY2更新應(yīng)力以及狀態(tài)變量，其中變量CTI（1）定義為彈性模量[E]，CTI（2）定義為泊松比[ν]。使彈性模量及泊松比作為線彈性材料的獨(dú)立常數(shù)，拉梅常數(shù)[λ]為

[λ=νE（1+ν）（1-2ν）.] （1）

另一組常用本構(gòu)常數(shù)體積模量[K]以及剪切模量[G]為

[K=E3（1-2ν），] （2）

[G=E2（1+ν）=μ.] （3）

通過廣義胡克定律可知

[σ=Cε，] （4）

式中：[σ]為柯西應(yīng)力；[C]為彈性剛度矩陣；[ε]為柯西應(yīng)變。式（2）和式（3）可表示為

[C=K+43GK-23GK-23G000K-23GK+43GK-23G000K-23GK-23GK+43G000 0 0 0G00 0 0 00G0  0 0 000G.] （5）

考慮到第[n]步以及第[n+1]步的時間間隔，即[tn， tn+1]，使用后退Euler法將本構(gòu)方程在時間間隔內(nèi)進(jìn)行離散，第[n+1]步的總應(yīng)變?yōu)?/p>

[εn+1=εn+Δεn+1]， （6）

式中：[Δεn+1]為第[n+1]步的應(yīng)變增量。

使用牛頓迭代法來更新第[n+1]步的應(yīng)力，即

[σn+1=σn+CΔεn+1]. （7）

3）KEY3更新剛度矩陣D，

[D=?Δσ?Δε]， （8）

式中：[Δσ]為柯西應(yīng)力增量；[Δε]為柯西應(yīng)變增量。

4）KEY4輸出結(jié)果及計(jì)算其他額外量。

3 求解器及子程序驗(yàn)證

3.1 桁架大變形結(jié)構(gòu)計(jì)算

為驗(yàn)證求解器的可用性以及求解精度，通過輸入模型文本文件計(jì)算了若干模型問題。圖2模型以及材料參數(shù)取自文獻(xiàn)[10]，為懸臂桁架結(jié)構(gòu)，對其進(jìn)行大變形分析，其中該結(jié)構(gòu)長寬比為9∶1，左端施加固定約束，右端中部節(jié)點(diǎn)施加豎直向下集中力載荷P，分別為50、75、100、125、150、175 N。桿件初始橫截面積A=10-4 m2，設(shè)置桿件彈性模量為2×109 N/m2。

通過設(shè)置不同載荷，并與商業(yè)軟件ABAQUS中打開幾何非線性開關(guān)進(jìn)行仿真模擬得出的數(shù)據(jù)進(jìn)行對比，得出大變形后載荷施加點(diǎn)橫向及縱向結(jié)果如圖3—圖4所示。從位移結(jié)果圖可以看出，利用求解器計(jì)算的結(jié)果與商業(yè)軟件ABAQUS計(jì)算結(jié)果最大誤差小于5%，符合工程上計(jì)算誤差，證明了求解器的準(zhǔn)確性。

3.2 子程序驗(yàn)證

為驗(yàn)證源代碼調(diào)用子程序進(jìn)行求解的可行性以及計(jì)算準(zhǔn)確性，本文編寫了基于廣義胡克定律的線彈性材料子程序，建立了8節(jié)點(diǎn)的三維實(shí)體單元模型，調(diào)用材料子程序來對其進(jìn)行單軸拉伸以及剪切狀態(tài)下的數(shù)值計(jì)算。模型如圖5所示，其中單軸拉伸時邊界條件為左側(cè)4節(jié)點(diǎn)固定，右側(cè)4節(jié)點(diǎn)沿Y軸方向拉伸0.01 m，模型彈性模量為2×1011 N/m2。

子程序計(jì)算的各點(diǎn)應(yīng)力結(jié)果如表7所示。通過表中數(shù)據(jù)可知，y方向應(yīng)力理論值應(yīng)為2×109 Pa，計(jì)算結(jié)果基本無誤差，其余方向應(yīng)力為計(jì)算容差且基本為0，計(jì)算容差正負(fù)符號符合邊界條件設(shè)置，證明了子程序的正確性以及使材料子程序嵌入求解器進(jìn)行求解的可行性。

4 結(jié)論

以ADINA84源代碼作為研究對象，對其編譯后求解器進(jìn)行大變形分析，通過子程序嵌入求解器進(jìn)行調(diào)用求解，將子程序計(jì)算結(jié)果與商業(yè)有限元軟件計(jì)算結(jié)果及理論值進(jìn)行對比，得出以下結(jié)論：

1）即使該源代碼誕生于40年前，但如今其依然能夠在保證精度的同時高效地完成線彈性、幾何非線性等有限元問題的求解；

2）該源代碼仍然可以通過用戶自定義材料模式完成二次開發(fā)工作，具備眾多商業(yè)有限元軟件核心功能；

3）基于ADINA84源代碼的編譯及子程序調(diào)用方法，可為國產(chǎn)有限元軟件開發(fā)提供借鑒和參考，這將為后續(xù)工作打下基礎(chǔ)。

參考文獻(xiàn)

[1] 貝斯.ADINA/ADINAT使用手冊：自動動態(tài)增量非線性分析有限元程序[M].趙興華，譯.北京：機(jī)械工業(yè)出版社，1986.

[2] 周世昌，劉從民，張錫清，等.利用結(jié)構(gòu)有限元程序ADINA、ADINAT進(jìn)行大型推力軸承的熱彈流計(jì)算[J].哈爾濱電工學(xué)院學(xué)報(bào)，1991（4）：351-358.

[3] 曾秋成.ADINA/ADINAT程序及其對車輪的應(yīng)力分析[J].鐵道車輛，1984，22（5）：29-38.

[4] 趙興華，徐福娣，梁醒培.ADINA/ADINAT程序述評[J].機(jī)械強(qiáng)度，1982，4（4）：15-20.

[5] 鄭結(jié).ADINA/ADINAT程序討論班[J].機(jī)械強(qiáng)度，1982，4（4）：72.

[6] 郭富強(qiáng)，李鋼鐵，霍軼珍，等.基于ADINA的河套灌區(qū)預(yù)制混凝土渠道保溫防凍脹效果數(shù)值模擬[J].干旱區(qū)資源與環(huán)境，2023，37（10）：84-92.

[7] 貢悅.基于數(shù)值模擬與實(shí)橋試驗(yàn)的拉索減振被動式阻尼器研究[D].大連：大連理工大學(xué)，2021.

[8] 謝建，權(quán)輝，田桂.基于ADINA的導(dǎo)彈起豎液壓缸有限元分析[J].起重運(yùn)輸機(jī)械，2011（7）：42-45.

[9] 陳志濤.卵形消化池流—固耦合系統(tǒng)的ADINA有限元分析[D].天津：天津大學(xué)，2007.

[10] 熊玉春，房營光.ADINA有限元軟件中材料本構(gòu)的二次開發(fā)[J].巖土力學(xué)，2008，29（8）：2221-2225，2240.

[11] 方瑞，陳岳坪.基于UG二次開發(fā)干涉避障仿真研究[J].廣西科技大學(xué)學(xué)報(bào)，2019，30（3）：64-69.

[12] 鄧?yán)誓?，黃曉霞，彭來，等.基于Revit二次開發(fā)的施工危險(xiǎn)源安全管理平臺研究與應(yīng)用[J].廣西科技大學(xué)學(xué)報(bào)，2018，29（1）：106-112.

[13] 徐武彬.AutoCAD二次開發(fā)技術(shù)在機(jī)械類課程教學(xué)中的應(yīng)用[J].廣西工學(xué)院學(xué)報(bào)，1999，10（1）：75-78.

[14] 張衛(wèi)東，岑威鈞.基于ADINA的黏土廣義塑性模型二次開發(fā)及驗(yàn)證[J].三峽大學(xué)學(xué)報(bào)（自然科學(xué)版），2016，38（6）：42-45.

[15] 汪仁和，李棟偉，王秀喜.改進(jìn)的西原模型及其在ADINA程序中的實(shí)現(xiàn)[J].巖土力學(xué)，2006，27（11）：1954-1958.

[16] 張宇鵬.拉壓不同模量材料的高性能算法開發(fā)及應(yīng)用[D].大連：大連理工大學(xué)，2016.

[17] 楊兵.ADINA幾何與材料雙重非線性源程序解析[D].南京：河海大學(xué)，2008.

Application and secondary development research of ADINA finite

element source code

QU Haoyu， MA Yingsong*， FAN Lizhong， GAO Sheng

（School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology，

Liuzhou 545616， China）

Abstract： As a large general finite element program， ADINA（automatic dynamic incremental nonlinear analysis）84 has a very powerful nonlinear computing ability and is widely used in various fields of engineering numerical simulation analysis. However， there are few studies on the analysis， application and secondary development of the source code. Through in-depth analysis of the source code， the material subroutine was written and the corresponding solver was compiled based on generalized Hooke's law， and the calculation accuracy of the solver and the material subroutine were compared and verified. The research provides a detailed description for the use of the ADINA source code and a reference for the secondary development of the constitutive model， and provides a new idea for the development of domestic special finite element software.

Keywords： ADINA source program; solver; geometric nonlinearity; secondary development

（責(zé)任編輯：于艷霞）

收稿日期：2023-10-18；修回日期：2023-11-07

基金項(xiàng)目：廣西科技基地和人才專項(xiàng)項(xiàng)目（桂科AD19110111）；廣西高校中青年教師科研基礎(chǔ)能力提升項(xiàng)目（2019KY0379）；中央引導(dǎo)地方科技發(fā)展資金項(xiàng)目（桂科ZY22096023）資助

第一作者：屈浩宇，在讀碩士研究生

*通信作者：馬迎松，博士，講師，研究方向：材料與結(jié)構(gòu)損傷破壞的宏細(xì)觀研究，E-mail：mys128@163.com