摘要

        上期我們介紹了LS-DYNA 新一代二次開發(fā)環(huán)境的編譯連接過程和新增功能，多個用戶子程序可以通過動態(tài)連接庫的方式同時動態(tài)加載。本文是其續(xù)篇，將以一個簡單的材料模型來演示在新的環(huán)境下的一個完整的開發(fā)過程，包括編譯，連接，動態(tài)加載，源程序跟蹤調試，以及模型驗證等環(huán)節(jié)。

引言

        LS-DYNA 作為一個大型的通用有限元程序，對于多重非線性的大規(guī)模問題的解決具有獨特的優(yōu)勢，在實際工程中也得到非常廣泛的應用。目前材料庫有300 種材料模型，其中多數都提供二維平面應力和三維應力兩個版本。LS-DYNA 提供完整的用戶材料模型的開發(fā)模板，讓用戶可以開發(fā)自己的材料模型。與一般的隱式算法相比，顯式有限元分析的時間步長很小，計算規(guī)模大，導致對用戶子程序的調用非常頻繁。LS-DYNA 為減少子程序的調用，內部采用批處理的方式調用用戶子程序，要求一次調用能處理幾百個單元，這也為用戶子程序實現矢量化計算提供了方便。因此，考慮到大變形，LS-DYNA 對用戶子程序的特殊要求也增加了用戶開發(fā)的復雜度。另外，對于一個對初次接觸LS-DYNA 的用戶來說，主程序的執(zhí)行碼不帶調試信息，較難在源程序上跟蹤調試，加大了二次開發(fā)中的程序查錯的難度。

        本文以一個簡單大變形下的線彈性材料模型為例，演示在新的開發(fā)環(huán)境下的完整的開發(fā)，調試和驗證過程。

線彈性材料物理模型

        1）應力應變關系

        在LS-DYNA 中應力和應變都是6 個分量，排列順序為

        其中E和v分別為楊氏模量和泊松比。

        2）驗證算例

        本文的驗證算例為只有一個8 節(jié)點體單元的模型，如圖一所示，長為L=2.0m，寬和高均為b=1.0m。加載條件為在X 方向單拉，而在Y 及Z 方向的位移為零。上述應力應變關系在小變形的情況下則簡化為

LS-DYNA 用戶子程序開發(fā)和調試

        1）UMAT 子程序的編譯和連接

        LS-DYNA 中的用戶材料號是從41 號到50 號，對應的第一級用戶入口子程序是dyn21.f 中的usrmat，受關鍵字*MAT_USER_DEFINED_MATERIAL_MODELS 控制。

        subroutine usrmat (lft,llt,cm,bqs,capa,eltype,mt,ipt,

       . npc,plc,crv,nnpcrv,rcoor,scoor,tcoor,nnm1,nip,ipt_thk)

        進入這個子程序后，再根據不同的單元類型選擇不同的第二級材料子程序

        ·urmathn: 體單元的三維材料模型

        ·urmats: 殼單元的二維平面應力材料模型

        ·urmatb, urmatd, urmatt: 三種不同的梁單元模型

        這三個不同子程序根據各自的單元特點對應力應變進行相應的第二級處理之后，再進入的第三級的用戶子程序umat41, umat42, … , umat50。這10 個子程序是標準的串行版本模板，演示不同類型的材料模型。一般的開發(fā)過程中，第一級和第二級入口程序都不需要改動,只需要從第三級這10 個子程序模板中選一個較為貼近的開始。本文選擇umat41 這個子程序。

        進入LS-DYNA 開發(fā)包的目錄后，進行如下幾個操作步驟：

        1. 把umat41 子程序從dyn21.f 中刪除，并復制到另一個新的源文件umat41.f

        subroutine umat41 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,

        1 temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject)

        2. 編輯開發(fā)包中的Makefile，把umat41.f的obj文件加到MY_OBJS變量中

        MY_OBJS = dyn21.o dyn21b.o init_dyn21.o couple2other_user.o dynrfn_user.o umat41.o

        3. 選擇編譯器

        原來的編譯器設置是和主程序一致的，LINUX系統(tǒng)一般是INTEL或者PGI的Fortran編譯器，這些商業(yè)編譯器的執(zhí)行代碼一般來說效率比較高。在LS-DYNA新的編譯環(huán)境下，用戶子程序的編譯器不要求和主程序一致，本文采用開源的gfortran來演示編譯過程。編譯環(huán)境為：

        Linux系統(tǒng)：OpenSUSE LEAP 42.1

        編譯器：gfortran 4.8.5

        MPI: platformmpi Community Edition 9.1.2

        ( http://www-03.ibm.com/systems/platformcomputing/products/mpi/ )

        將Makefile中的編譯變量設置為

MY_FLAG = -g -fPIC -fcray-pointer -I/opt/platform_mpi/include

FC = /usr/bin/gfortran

LD = /usr/bin/gfortran -shared

export MPI_F77 := /usr/bin/gfortran

MY_TARGET = gnu.so

        其中-g是讓編譯的用戶模塊帶源程序的調試跟蹤信息。這些變量的詳細解釋請參閱上期的“LS-DYNA用戶子程序的編譯和連接”一節(jié)。

        4.用make命令編譯，生成gnu.so，就完成了編輯和連接。

        2）UMAT子程序的調用

        上面編譯好的gnu.so可以做為開發(fā)好的用戶模塊配合模型使用。這個模塊和LS-DYNA主執(zhí)行程序是分開的，即使將來LS-DYNA主程序的版本升級也不影響這個模塊。調用的方法是在模型的.k文件里面加入三行

*MODULE_LOAD

myumat41

gnu.so

其中：第一行是關鍵字，第二行是這個模塊在這個模型的id，第三行是這個模塊的編譯后文件。然后就可以按照原來的方法執(zhí)行LSDYNA主程序就可以了。這個關鍵字有很多匹配規(guī)則，詳見上期的“LS-DYNA用戶子程序的動態(tài)連接庫的調用”一節(jié)。本文演示的是MPP版本的主程序，單個單元模型只能用一個CPU來運行：

mppdyna i=demo.k

        3）UMAT子程序的跟蹤調試

        當子程序運行遇到問題的時候，簡單直接方法是用打印命令，與LS-DYNA的59號文件對應的是message文件，對于MPP程序，每個CPU都有一個自己的message文件，因此打印方法不容易混亂，也很方便。比如：

WRITE(59,*)’sig=’,sig(1),sig(2),sig(3)

WRITE(59,*)’hsv=’,hsv(1),hsv(2)

        有些情況下，還是要進入到源程序里面，在源程序上進行跟蹤調試。本文以gdb為例，啟動調試程序，進行以下步驟：

        1. 調入主程序

gdb mppdyna

        2. 設置斷點

        b umat41

        (此時會顯示umat41 不存在，可能要用set breakpoint pending on 激活在調用時補設斷點)

        3. 運行程序

        r i=demo.k

        4. 程序在進入umat41 后，就會停下來。比如，打印變量

        p cm(1)

        這是* MAT_USER_DEFINED_MATERIAL_MODELS 卡片輸入的第一個材料常數P1。

材料模型的驗證

        根據前面定義的物理模型，利用LS-Prepost 建立一個有限元模型，包含一個八節(jié)點的實體單元，長度L=2m，寬和高為b=1m。并用LS-Prepost 施加所有的邊界條件和給定位移以及材料屬性后，有限元模型如圖二所示：

        所有的節(jié)點都固定，只有2，4，6，8 節(jié)點在X 方向有指定速度為1.0m/s。分析時間為1s， 則節(jié)點位移和工程應變時間曲線為：

        與圖三的結果吻合得很好。

        材料模型為線彈性，E=150x10⁹ 及v=0.25，則應力時間曲線為：

        σ(t)=180 x10⁹ε(t)=90x10⁹t

        而圖四中計算結果給出的最大應力則只有σmax =73 x10⁹Pa  

        原因是LS-DYNA 中的應變都是真實應變，而不是上面計算的工程應變。真實應變的計算方法為：

        得到最大真實應變e_max=ln(1.5)=0.405

        代入應力計算公式可以得到最大真實應力的理論值是73 x10⁹Pa，與有限元的分析結果完全吻合。同時計算結果也表明，當時間t=0.001, 工程應變ε(0.001)=0.005，僅為千分之五， 而此時的真實應力為σ(0.001)=89.7x10⁹t, 與線性公式相差0.3%，而當時間t=0.2時，工程應變達到百分之十，ε(0.2)=0.1，真實應力為σ(0.2)=171 x10⁹t, 與線性公式相差4.7%。

        在開發(fā)LS-DYNA 的用戶子程序的時候，一定要考慮大變形情況下的本構關系。一般來說， 把線性小變形的本構關系直接放進LS-DYNA，真實應力會有很大的差別。這點與一般的小變形下的本構模型開發(fā)有很大的不同，也加大了LS-DYNA 二次開發(fā)的難度。

作者簡介

*韓志東/Zhidong Han博士1998年畢業(yè)于清華大學計算固體力學專業(yè)，于2011年加入LSTC。他目前從事材料損傷斷裂分析及厚殼單元等方面研發(fā)。