1 uel簡要介紹

1.1 uel變量

在uel中，主要需要更新amatrix（雅克比矩陣）和rhs（右端殘值項(xiàng)），如果有狀態(tài)變量的話，也需要更新狀態(tài)變量svars.

1.2 inp文件修改

并且要成功使用uel的話，需要對(duì)inp文件做以下修改，一一說明。

（1） 需要首先定義uel，加入以下keyword以及data line(s)

各變量的含義如下：

該keyword下的data line(s)表示uel激活的自由度。

（2） 定義完uel之后，可以使用該uel，keyword以及data line(s)如下：

（3） 接著需要給uel賦予性質(zhì)

以C3D8單元計(jì)算一個(gè)各向同性材料為例，uel需要的性質(zhì)為材料的楊氏模量和泊松比。

1.3 uel結(jié)果的可視化

由于Abaqus不支持uel結(jié)果的可視化，因此我們可以利用umat來輔助進(jìn)行可視化。使用Abaqus的標(biāo)準(zhǔn)單元覆蓋uel，即標(biāo)準(zhǔn)單元與uel共享節(jié)點(diǎn)，但是材料使用umat（設(shè)置一個(gè)很小的剛度），關(guān)鍵步驟在于使用common block將uel中的svars數(shù)據(jù)傳遞到umat中的statev中。

2 有限元離散方程

2.1 形函數(shù)

本部分記錄八節(jié)點(diǎn)實(shí)體單元（C3D8）的剛度矩陣和內(nèi)力列陣等 ，熟悉此部分可跳過。

該單元的節(jié)點(diǎn)編號(hào)順序示意圖如下：

單元的形函數(shù)矩陣為：

參數(shù)坐標(biāo)系下的形函數(shù)為：

因此形函數(shù)對(duì)參數(shù)坐標(biāo)的導(dǎo)數(shù)為：

2.2 雅克比矩陣

利用形函數(shù)對(duì)參數(shù)坐標(biāo)的導(dǎo)數(shù)矩陣和節(jié)點(diǎn)的物理坐標(biāo)信息可以求得雅克比矩陣：

具體公式為：

2.3 B矩陣

利用雅克比矩陣可以求得形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)：

從而可以組裝B矩陣為：

2.4 單元?jiǎng)偠染仃嚭蛢?nèi)力列陣

單元的剛度矩陣和內(nèi)力公式為：

uel的主要目的就是更新這兩項(xiàng)，其中和的計(jì)算需要用到具體的本構(gòu)方程。在此使用各向同性線彈性本構(gòu)。

2.5 高斯積分

上述積分的計(jì)算需要使用數(shù)值積分，一維高斯積分點(diǎn)的坐標(biāo)以及權(quán)系數(shù)如下：

采用2x2x2的高斯積分點(diǎn)。

3 B bar 算法

B-bar算法將應(yīng)變修正為:

式中:

體應(yīng)變插值B矩陣為:

記為:

也有:

因此將B矩陣修正為:

4 測(cè)試

4.1 一個(gè)單元測(cè)試

對(duì)一個(gè)單元進(jìn)行單軸拉伸測(cè)試，借助umat進(jìn)行可視化時(shí)，在代碼中設(shè)置參數(shù)為：

需要根據(jù)job的單元總數(shù)手動(dòng)修改該參數(shù)。

位移結(jié)果為：

應(yīng)變結(jié)果為：

應(yīng)力結(jié)果為：

4.2 帶孔板單軸拉伸測(cè)試

位移結(jié)果對(duì)比為（左列為Abaqus計(jì)算結(jié)果，右列為uel計(jì)算結(jié)果）：

應(yīng)變的對(duì)比結(jié)果為（左列為Abaqus計(jì)算結(jié)果，右列為uel計(jì)算結(jié)果）：

應(yīng)力的對(duì)比結(jié)果為（左列為Abaqus計(jì)算結(jié)果，右列為uel計(jì)算結(jié)果）：

uel源代碼

2024/03/27更新:已加入B-bar算法進(jìn)行剛度矩陣的修正,可看到uel的計(jì)算結(jié)果已于Abaqus計(jì)算結(jié)果一致