大家好，第一次在技術(shù)鄰平臺發(fā)貼，如有不當(dāng)歡迎指出。

因為課題原因，需要用到cohesive單元連接兩殼單元的邊界，然而ABAQUS中的cohesive單元只有用于連接實體（solid）或者平面單元（平面應(yīng)力或平面應(yīng)變）的coh單元，在導(dǎo)師的催促下就只能從零開始自學(xué)自定義單元的內(nèi)容，最后也終于基本實現(xiàn)這個目標(biāo)。自學(xué)的過程很辛苦，其中很多東西考慮的也并不是很充分，理解的也不一定到位，但結(jié)果實現(xiàn)了目標(biāo)，還是很開心的。

當(dāng)然這個過程中很感謝snowwave02和借風(fēng)一尺兩位的幫助，學(xué)術(shù)方面大家多交流才會進(jìn)步。

話不多說，進(jìn)入正題。其實這種單元在DYNA中是有的，但是因為課題原因希望進(jìn)行對coh單元本構(gòu)更多的二次開發(fā)，且DYNA的二次開發(fā)我不是很熟悉，因此我就直接把DYNA的理論搬過來用了。但是鑒于二者還是有所不同，因此我在一些細(xì)節(jié)部分進(jìn)行了個人的修改。這次的介紹主要理論源自于一篇DYNA會議上的文章：Edge-to-edge Cohesive Shell Elements in LS-DYNA, 作者是Jesper Karlsson, 如果有興趣可以找原文查看。

單元的基本形式如上圖所示，四個節(jié)點，N1,N2用于和上邊的殼單元（ABAQUS關(guān)鍵字S4）相連，N3，N4則和下邊的殼連接。殼的厚度為t，采用常規(guī)的四個積分點（圖中畫“×”的位置），積分點的等參坐標(biāo)也是非常常見的： (±1/√3,±1/√3)。對于每個節(jié)點均開放六個自由度，即三個平動三個轉(zhuǎn)動。

接下來是運動學(xué)的定義：

d代表分離量，即traction-separation法則中的分離距離，xt和xb分別代表上下邊在對應(yīng)等參坐標(biāo)下的坐標(biāo)位置，xi（i=1，2，3，4）代表節(jié)點的坐標(biāo)，η 和  ξ 是等參坐標(biāo)。按照這種方式算出來在積分點的分離量d如下圖紅線所示意，因為有四個積分點，恰好有四段d，這也是很清晰的。R則是旋轉(zhuǎn)矩陣，因為traction-separation需要在局部坐標(biāo)系下討論，因此需要旋轉(zhuǎn)矩陣溝通局部與全局坐標(biāo)系。ni（i=1，2，3，4）代表對應(yīng)節(jié)點處的法向方向。

本構(gòu)關(guān)系方面，就比較簡單：

其中d是剛剛定義的分離量，T則是cohesive法則中的牽引力，K是剛度，D是損傷系數(shù)，因此只需要完成D的更新，就可以實現(xiàn)本構(gòu)。不同的traction-separation法則D的更新有所不同，這篇帖子主要不是講本構(gòu)的更新，因此這一部分就略過。

VUEL中需要用戶代碼定義三個內(nèi)容：質(zhì)量矩陣、節(jié)點內(nèi)力以及臨界穩(wěn)定步長，首先討論節(jié)點力。

定義節(jié)點力一般需要從虛功率方程出發(fā)，首先我們寫出虛功率方程：

上式中，d是之前介紹的分離量，T為牽引力，f是節(jié)點內(nèi)力，m是節(jié)點內(nèi)轉(zhuǎn)矩，x是節(jié)點坐標(biāo)，omiga是角速度。

另一方面，上面虛功率方程的左端可以展開，直接帶入帖子最前面運動學(xué)的定義：

上式中忽略最后一項，因為認(rèn)為旋轉(zhuǎn)矩陣的變化是很小的，將xt和xb分別對時間求導(dǎo)（二者的表達(dá)式上文已經(jīng)給出）再帶入，對比虛功率方程右端以及剛剛展開式子的右端，就可以得到節(jié)點內(nèi)力以及節(jié)點彎矩的表達(dá)式：

（表達(dá)式太多就不打全了，可以自行推導(dǎo)）

當(dāng)然，要找到上面積分的原函數(shù)是很困難的，帖子最開始說的數(shù)值積分就是需要用在這個時候的！由高斯積分法，函數(shù)在某個區(qū)域的積分可以看作在區(qū)域內(nèi)某幾個點對應(yīng)函數(shù)值的求和，寫成表達(dá)式的話有：

其中

有了數(shù)值積分的工具，我們就可以把上面求得的積分式算出來，我把他們?nèi)苛谐鰜恚?/p>

觀察上面的這些式子，我們看一下哪些值我們還沒有。首先，T是通過本構(gòu)關(guān)系求得的，已經(jīng)更新完畢；積分點的等參坐標(biāo)都是給定的；t是參數(shù)，也是知道的；所以我們沒有的量是旋轉(zhuǎn)矩陣R和法向量n。

首先講講節(jié)點的法向量n，在下圖中用紅色箭頭示意，這里面n默認(rèn)是單位向量。

ABAQUS在更新是都會給出每個節(jié)點的轉(zhuǎn)角增量，那么我們從這個轉(zhuǎn)角增量，根據(jù)羅德里格旋轉(zhuǎn)公式，就可以得到新的法向量：

這樣我們每一步都會更新一個n，然后我們把他存到歷史變量中，它就會隨著有限元計算一步一步更新了

好了，下面介紹旋轉(zhuǎn)矩陣R的更新：

這樣我們就有了R：

同樣把R存到狀態(tài)變量中，一步步更新

至此我們就有了f和m中所有需要的量，也就可以完成節(jié)點內(nèi)力（彎矩）的所有更新

然后是質(zhì)量矩陣，我的質(zhì)量矩陣就比較隨意，采用的是集中質(zhì)量矩陣，意思是只有對角線上有非零值，由于質(zhì)量守恒，因此質(zhì)量矩陣應(yīng)該是個常量，所以我就在對角線上隨意給了一個量（意思是對角線上都是這同一個值，其他地方都是0），沒想到最后也能跑通，關(guān)于質(zhì)量矩陣怎么給，也想和大家更多交流。

最后是臨界時間步長，我也是參考有限元書上的公式自己定義了一個量：

其中beta是系數(shù)，按照書中內(nèi)容我取了0.95，E和rho分別代表模量和密度。

以上就將所有量更新完畢，這個單元就可以用了。

展示一下效果：

兩個大的單元都是純彈性的殼單元（S4），打“×”的地方就是剛剛寫的VUEL，還是實現(xiàn)了預(yù)期的功能了的，當(dāng)然后續(xù)還有更多的調(diào)試，也歡迎大家多交流