復(fù)合材料機(jī)械連接結(jié)構(gòu)分析中的接觸算法適用性研究

李沛城，趙美英，侯赤

(西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)

0 引 言

現(xiàn)階段復(fù)合材料機(jī)械連接結(jié)構(gòu)仍然大量存在于飛機(jī)結(jié)構(gòu)之中，而該類結(jié)構(gòu)受載時(shí)存在典型的非線性接觸效應(yīng)，早期的學(xué)者通過解析法對(duì)此展開了研究，取得了一定成果。F.K.Chang等[1]和H.T.Hahn等[2]分別采用孔邊余弦分布載荷假設(shè)和剪切非線性本構(gòu)關(guān)系研究了含孔層合板銷釘加載條件下的應(yīng)力分布情況。然而，銷釘與層合板接觸邊界上的接觸內(nèi)力受銷釘彈性、層合板鋪層方式以及摩擦等多方面因素影響，所以余弦分布載荷假設(shè)是存在缺陷的。陳浩然等[3]提出了一種處理非線性接觸問題的增量有限元混合法，并分析了銷釘彈性、層合板鋪層以及摩擦效應(yīng)對(duì)復(fù)合材料機(jī)械連接結(jié)構(gòu)孔邊非線性接觸應(yīng)力分布的影響。

目前，隨著商用有限元軟件的快速發(fā)展，它們被大量運(yùn)用在含接觸工程結(jié)構(gòu)的應(yīng)力分析工作中。張寶生等[4]采用Marc軟件對(duì)圓柱滾子接觸問題進(jìn)行求解，發(fā)現(xiàn)接觸區(qū)域內(nèi)網(wǎng)格數(shù)量越少，計(jì)算結(jié)果誤差越大。龔思楚等[5]認(rèn)為采用固定約束模擬典型連接件的緊固件使得結(jié)構(gòu)局部剛度過硬，容易導(dǎo)致連接部位計(jì)算結(jié)果失真；而當(dāng)采用Abaqus軟件建立接觸對(duì)模擬連接緊固件后，所得計(jì)算結(jié)果與試驗(yàn)結(jié)果較為吻合。陳向明等[6]基于Abaqus自定義了一種只能傳遞壓力而不能傳遞拉力的粘接元，將其設(shè)置在兩接觸對(duì)之間來傳遞接觸壓力，有效地解決了復(fù)合材料接頭強(qiáng)度分析中由接觸帶來的收斂困難問題。

對(duì)于復(fù)合材料機(jī)械連接結(jié)構(gòu)孔邊應(yīng)力分布及連接強(qiáng)度的研究當(dāng)前大多借助于數(shù)值分析方法進(jìn)行[6-9]，但由于材料以及機(jī)械連接結(jié)構(gòu)的復(fù)雜性使得計(jì)算結(jié)果影響因素多且計(jì)算分散度大,其中接觸算法的選取對(duì)有限元計(jì)算結(jié)果會(huì)產(chǎn)生明顯影響。為了模擬連接結(jié)構(gòu)中連接件與被連接件之間的載荷傳遞, 需要在二者接觸面上設(shè)置接觸邊界，常用的接觸邊界建立方法分為兩類：一是在接觸邊(面)上建立接觸單元，例如Nastran中的Gap元、Slidline單元等；二是在接觸邊(面)上定義非線性接觸邊界條件，例如Marc中的柔性體接觸，Abaqus中的點(diǎn)面、面面接觸等。無(wú)論以上述何種方式設(shè)置接觸，在有限元計(jì)算中都會(huì)依據(jù)特定的接觸算法進(jìn)行接觸狀態(tài)檢測(cè)和施加接觸約束，且對(duì)于不同的分析模型，這些接觸算法的適用性存在差異。

本文針對(duì)復(fù)合材料機(jī)械連接結(jié)構(gòu)中螺栓與層合板間的接觸特點(diǎn)，對(duì)一些典型的接觸算法在復(fù)合材料機(jī)械連接應(yīng)力分析中的準(zhǔn)確性和適用性進(jìn)行對(duì)比研究，闡述其原理與影響因素，以期為此類結(jié)構(gòu)的有限元分析工作提供參考。

1 接觸分析基本流程

含接觸的有限元分析是一個(gè)依賴于時(shí)間并伴隨結(jié)構(gòu)邊界非線性演化的過程，因此接觸問題一般采用增量法求解。由于接觸面的范圍和接觸狀態(tài)事先未知，在求解過程中需要進(jìn)行多次迭代來尋求結(jié)構(gòu)的平衡狀態(tài)。具體而言包括三個(gè)主要步驟：(1) 根據(jù)前一增量步的結(jié)果和當(dāng)前增量步給定的載荷條件，通過接觸狀態(tài)檢測(cè)，設(shè)定此增量步初次迭代步的接觸區(qū)域和接觸狀態(tài)；(2) 對(duì)假定發(fā)生接觸的點(diǎn)施加接觸力或接觸位移約束，并引入總體協(xié)調(diào)方程進(jìn)行求解；(3) 利用計(jì)算結(jié)果對(duì)假定的接觸狀態(tài)進(jìn)行檢查。如果接觸面上每一點(diǎn)都不違反假定狀態(tài)，則完成本增量步的求解并轉(zhuǎn)入下一增量步的求解：否則修改接觸狀態(tài)，回到步驟(2)進(jìn)行下一次迭代求解。

以上步驟的核心技術(shù)包括接觸檢測(cè)技術(shù)和接觸后節(jié)點(diǎn)約束的施加技術(shù)。前者用于確定發(fā)生接觸的區(qū)域；后者用于對(duì)發(fā)生接觸的點(diǎn)在系統(tǒng)方程中引入所需要滿足的位移和力邊界條件。

2 接觸檢測(cè)技術(shù)

接觸檢測(cè)技術(shù)用于判定發(fā)生接觸的單元和節(jié)點(diǎn)、發(fā)生接觸的方向和局部坐標(biāo)系以及接觸發(fā)生的判據(jù)；根據(jù)接觸對(duì)的判定過程是在計(jì)算前還是在計(jì)算中進(jìn)行，接觸檢測(cè)技術(shù)可以分為小滑移接觸檢測(cè)技術(shù)和有限滑移接觸檢測(cè)技術(shù)。

2.1 小滑移接觸檢測(cè)技術(shù)

小滑移接觸中接觸對(duì)的判定是在計(jì)算前進(jìn)行，計(jì)算中接觸對(duì)將不再變化。當(dāng)有接觸單元存在時(shí)，接觸對(duì)即為接觸單元的相應(yīng)節(jié)點(diǎn)，例如Nastran中的Gap元(如圖1所示[10])；在沒有明確的接觸單元存在時(shí)，接觸對(duì)是由節(jié)點(diǎn)或其投影點(diǎn)與單元邊界在某一局部坐標(biāo)系中的相對(duì)位置確定。

圖1 Gap元及其局部坐標(biāo)系

圖1所示為Nastran和Marc提供的用于模擬點(diǎn)點(diǎn)接觸的Gap元及其局部坐標(biāo)系[10]，Gap元指定的接觸對(duì)是Gap元的GA和GB兩節(jié)點(diǎn)，當(dāng)Gap元被定義好后，相應(yīng)的局部坐標(biāo)系也將固定不變，因此接觸對(duì)和沿著坐標(biāo)系x軸的接觸方向在計(jì)算過程中均保持不變。判定接觸發(fā)生后，Nastran中Gap元根據(jù)節(jié)點(diǎn)的相對(duì)位移量對(duì)接觸節(jié)點(diǎn)施加接觸力；而Marc中的Gap元?jiǎng)t采用拉格朗日乘子法來施加接觸約束。

除了利用接觸單元，還可以根據(jù)節(jié)點(diǎn)與單元邊界之間的相對(duì)位置判定接觸對(duì)，一種典型的從面節(jié)點(diǎn)與主面接觸關(guān)系判定的示意圖如圖2所示[11]。在計(jì)算開始之前的模型分析過程中，首先確定主面單元兩端節(jié)點(diǎn)的法向量，由單元節(jié)點(diǎn)的法向量與單元的插值函數(shù)就可以共同確定出單元中每個(gè)位置的法向量。對(duì)于與主面接觸的從節(jié)點(diǎn)，可以在主面單元中找到一個(gè)對(duì)應(yīng)的投影點(diǎn)，也叫錨點(diǎn)，主面單元在錨點(diǎn)處的法向量恰好通過從節(jié)點(diǎn)，再根據(jù)錨點(diǎn)和錨點(diǎn)處的法線方向，便可以確定與從節(jié)點(diǎn)對(duì)應(yīng)的切平面。由從節(jié)點(diǎn)與切平面的相對(duì)位置即可判斷接觸是否發(fā)生。

圖2 小滑移接觸檢測(cè)示意圖

在小滑移接觸檢測(cè)中通常只考慮切平面的轉(zhuǎn)動(dòng)和默認(rèn)錨點(diǎn)的位置不改變，且一旦確立切平面后，從節(jié)點(diǎn)和切平面的接觸對(duì)應(yīng)關(guān)系也不再發(fā)生變化。

2.2 有限滑移接觸檢測(cè)技術(shù)

有限滑移中接觸對(duì)的判定是在計(jì)算中不斷檢測(cè)并更新的。因此接觸對(duì)關(guān)系會(huì)發(fā)生變化，當(dāng)接觸面上有切向滑移位移存在時(shí)，采用此種檢測(cè)技術(shù)與實(shí)際情況更為吻合。

Marc在處理柔性體接觸問題時(shí)所提供的接觸狀態(tài)判定方法如圖3所示[12]，其基本思想基于接觸邊(面)的接觸容限，特點(diǎn)在于它并不以被檢測(cè)節(jié)點(diǎn)和接觸段上某個(gè)節(jié)點(diǎn)在某個(gè)特定方向上的距離作為判斷依據(jù)，而是在接觸段兩側(cè)設(shè)定一個(gè)接觸容限區(qū)域，只要節(jié)點(diǎn)落入這個(gè)區(qū)域即認(rèn)為節(jié)點(diǎn)與對(duì)應(yīng)的接觸段發(fā)生接觸。在使用牛頓-拉普森算法時(shí)，Marc將在每個(gè)增量步的迭代中進(jìn)行穿透檢測(cè)。

圖3 接觸距離容限示意圖

此外，圖2所示接觸檢測(cè)方式如果設(shè)定為在計(jì)算中由程序跟蹤從節(jié)點(diǎn)的位置，并不斷更新投影點(diǎn)和切平面，從而確定新的接觸狀態(tài)，此時(shí)也就轉(zhuǎn)變?yōu)橛邢藁平佑|檢測(cè)。

3 接觸約束的施加技術(shù)

常用的接觸約束施加方法包括拉格朗日乘子法、罰函數(shù)法、直接約束法[13-15]。

3.1 拉格朗日乘子法

拉格朗日乘子法是通過拉格朗日乘子施加接觸體必須滿足的非穿透約束條件，是一種帶約束極值問題的描述方法，它是把約束條件加在一個(gè)系統(tǒng)中最完美的數(shù)學(xué)描述。在迭代求解中，應(yīng)用拉格朗日乘子法時(shí)的系統(tǒng)控制方程一般形式為

(1)

該方法雖然增加了系統(tǒng)變量數(shù)目，且需要在數(shù)值方案中處理非正定系統(tǒng)，需要額外的操作才能保證計(jì)算精度，但計(jì)算中無(wú)需用戶指定穿透剛度，可以避免由于穿透剛度選擇不當(dāng)而帶來的計(jì)算不收斂問題。

3.2 罰函數(shù)法

罰函數(shù)法的原理是一旦接觸區(qū)域發(fā)生穿透，罰函數(shù)便夸大這種誤差的影響，從而使系統(tǒng)的求解(滿足力的平衡和位移的協(xié)調(diào))無(wú)法正常實(shí)現(xiàn)。從物理意義上講，用罰函數(shù)法施加接觸約束相當(dāng)于在物體之間指定彈性變形的區(qū)域，其作用類似非線性彈簧[16-17]，該方法不增加未知量數(shù)目，但會(huì)增加系統(tǒng)矩陣帶寬。在實(shí)際的計(jì)算中，罰函數(shù)通常不是一個(gè)恒定的常量， Abaqus常用的指數(shù)型罰函數(shù)接觸約束關(guān)系如圖4所示[11]。

圖4 指數(shù)型罰函數(shù)與穿透量關(guān)系

3.3 直接約束法

用直接約束法處理接觸問題的原理是追蹤物體的運(yùn)動(dòng)軌跡，—旦探測(cè)出發(fā)生接觸，便將接觸所需的運(yùn)動(dòng)約束和節(jié)點(diǎn)力作為邊界條件直接施加在產(chǎn)生接觸的節(jié)點(diǎn)上。例如在分析柔性體與剛體接觸的時(shí)候，約束條件是不允許有穿透發(fā)生，此時(shí)系統(tǒng)直接給柔性體中有接觸的節(jié)點(diǎn)施加對(duì)應(yīng)的位移邊界條件，使得節(jié)點(diǎn)的位移增量和與之發(fā)生接觸的剛性面變化一致，每一迭代步的位移約束量滿足如下方程式：

(2)

式中：n1,n0為不同迭代步的剛體表面法線方向向量；v為預(yù)估的剛體表面位移。

該方法不需要增加特殊的界面單元，也不增加系統(tǒng)自由度數(shù)，但由于接觸關(guān)系的變化會(huì)增加系統(tǒng)矩陣帶寬。另外，直接約束法的關(guān)鍵在于準(zhǔn)確預(yù)測(cè)施加在從節(jié)點(diǎn)的約束量。對(duì)于柔性體與剛體接觸的情況，在直接約束法中由于節(jié)點(diǎn)位移約束是根據(jù)接觸面的預(yù)知位移或速度以及表面形狀來確定的，因此，當(dāng)接觸面網(wǎng)格越細(xì)致，表面形狀就越準(zhǔn)確，相應(yīng)的位移約束也就越準(zhǔn)確，計(jì)算精度也會(huì)更高。

4 接觸算法適用性對(duì)比算例

為了比較不同接觸算法在復(fù)合材料機(jī)械連接結(jié)構(gòu)應(yīng)力分析中的適用性，以Wang W C等[18]開展的單向板銷釘連接應(yīng)力試驗(yàn)中的試驗(yàn)件為計(jì)算對(duì)象，利用有限元方法對(duì)單向板銷釘連接試驗(yàn)件進(jìn)行孔邊應(yīng)力分布數(shù)值計(jì)算(模型參數(shù)詳見文獻(xiàn)[18])，并將計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)[18]進(jìn)行對(duì)比。

4.1 接觸檢測(cè)技術(shù)適用性對(duì)比算例

在接觸檢測(cè)技術(shù)中，由Gap元所定義的點(diǎn)點(diǎn)接觸是用戶在建模時(shí)就指定將要發(fā)生接觸的節(jié)點(diǎn)對(duì)，并依據(jù)單元兩節(jié)點(diǎn)的相對(duì)位置判定接觸狀態(tài)，屬于小滑移接觸檢測(cè)，不需要進(jìn)行接觸探測(cè)，計(jì)算規(guī)模較小；而Marc中基于接觸容限的接觸檢測(cè)可以設(shè)定為在計(jì)算中實(shí)施接觸檢測(cè)和接觸對(duì)更新，因此屬于有限滑移接觸檢測(cè)。為了對(duì)比二者的適用性，分別建立相應(yīng)的有限元模型計(jì)算釘孔邊緣徑向應(yīng)力并與試驗(yàn)值對(duì)比，如圖5所示。

圖5 接觸檢測(cè)對(duì)孔邊徑向應(yīng)力分布的影響

由于試驗(yàn)件沿中心線方向施加拉伸載荷，在以孔中心為原點(diǎn)的圓柱坐標(biāo)系中，隨著孔邊各點(diǎn)與中心線夾角由90°向0°趨近，它們對(duì)應(yīng)的徑向應(yīng)力總體上呈逐漸增大的趨勢(shì)。但值得注意的是，孔邊最大徑向應(yīng)力并沒有出現(xiàn)在0°位置，而是在10°位置附近。M.W.Hyer等[19]對(duì)這種現(xiàn)象進(jìn)行了研究并認(rèn)為這是由接觸中的粘著效應(yīng)和摩擦力共同作用所導(dǎo)致，而在有限元計(jì)算中，因?yàn)闊o(wú)法模擬接觸的粘著效應(yīng)，所以接觸面之間必定存在一定程度的切向滑移量。

圖5中采用基于Gap元的小滑移檢測(cè)技術(shù)所得計(jì)算結(jié)果準(zhǔn)確度低于有限滑移接觸模型所得結(jié)果，其原因就是連接件與釘孔的接觸面上有切向滑移存在，而Gap元的接觸對(duì)不會(huì)隨著接觸面滑移而改變，所以此時(shí)再沿著原Gap元的方向判斷接觸狀態(tài)已失效。

為了進(jìn)一步說明Gap元的計(jì)算精度與切向位移量有關(guān)，圖5分別給出了1 mm和2 mm位移載荷作用下釘孔周圍的徑向應(yīng)力分布，由于位移載荷大小會(huì)影響釘孔的變形量，導(dǎo)致接觸面上出現(xiàn)不同的切向位移，所以在不同位移載荷下計(jì)算所得徑向應(yīng)力比也不相同。

由此可知，對(duì)于螺栓連接結(jié)構(gòu)應(yīng)盡量采用有限滑移接觸檢測(cè)的模型。

4.2 接觸約束施加技術(shù)適用性對(duì)比算例

采用拉格朗日乘子法施加接觸約束通常用于接觸狀態(tài)可以預(yù)知的情況之下，此時(shí)所得解非常精確。但對(duì)于復(fù)合材料螺栓連接結(jié)構(gòu)，由于切向滑移的存在，接觸檢測(cè)的誤差將會(huì)完全引入分析過程，嚴(yán)重影響計(jì)算結(jié)果，因此下文只針對(duì)罰函數(shù)法和直接約束法進(jìn)行對(duì)比分析。

罰函數(shù)法將在模型中引入較大的穿透剛度，且罰函數(shù)只有在穿透發(fā)生后才生效，因此會(huì)出現(xiàn)節(jié)點(diǎn)穿透的情況；而直接約束法則是將接觸的位移約束轉(zhuǎn)化到接觸段的局部坐標(biāo)系中，作為邊界條件直接施加到節(jié)點(diǎn)上，使用直接約束法不會(huì)產(chǎn)生節(jié)點(diǎn)的穿透現(xiàn)象。本文分別采用罰函數(shù)法和直接約束法，以相同的有限元模型計(jì)算釘孔邊的徑向應(yīng)力，計(jì)算結(jié)果如圖6所示。

圖6 接觸約束對(duì)孔邊徑向應(yīng)力分布的影響

從圖6可以看出：在孔邊30°范圍以內(nèi)二者計(jì)算結(jié)果幾乎一致，但是在約30°的位置由直接約束法所得計(jì)算結(jié)果出現(xiàn)了明顯的奇異點(diǎn)，從而導(dǎo)致30°至40°范圍內(nèi)的計(jì)算結(jié)果與罰函數(shù)法差異較大。出現(xiàn)奇異點(diǎn)的原因在于直接約束法的求解精度依賴于對(duì)發(fā)生接觸的節(jié)點(diǎn)所應(yīng)滿足的位移條件的預(yù)估準(zhǔn)確度。在本例中，此位移條件由銷釘?shù)奈灰坪徒佑|表面形狀所決定，當(dāng)銷釘固定時(shí)，其位移保持為0。而接觸表面形狀則是利用樣條函數(shù)和單元形函數(shù)共同描述的，因此與實(shí)際幾何出現(xiàn)差異，如果增加接觸面的單元密度，就能改善計(jì)算精度。

保持圖6對(duì)應(yīng)的有限元模型其余設(shè)置不變，僅將網(wǎng)格密度增加一倍后分別采用直接約束法和罰函數(shù)法所得計(jì)算結(jié)果如圖7所示，證明了網(wǎng)格適當(dāng)加密后，局部的奇異點(diǎn)得到有效消除，計(jì)算精度進(jìn)一步提高。

圖7 接觸約束對(duì)孔邊徑向應(yīng)力分布的影響(網(wǎng)格加密)

由此可見，雖然直接約束法可以不需要用戶輸入接觸剛度等數(shù)據(jù)就能獲得收斂的結(jié)果，但是需要用戶建立與實(shí)際結(jié)構(gòu)嚴(yán)格一致的模型和精度足夠的網(wǎng)格，導(dǎo)致有限元建模的前處理工作相對(duì)繁瑣；罰函數(shù)法對(duì)模型精度的要求則沒有直接約束高，只需要用戶選擇合適的罰函數(shù)就可以保證結(jié)果收斂。

5 結(jié) 論

(1) 復(fù)合材料機(jī)械連接結(jié)構(gòu)中連接件與釘孔之間存在著一定的切向滑移，因此盡量采用有限滑移接觸檢測(cè)方式，避免使用Gap元。

(2) 接觸約束施加方法中，拉格朗日乘子法雖是施加約束最為準(zhǔn)確的一種方法，但是需要準(zhǔn)確判定接觸狀態(tài)，只適用于無(wú)滑移時(shí)的接觸狀況，因此不適用于易產(chǎn)生滑移的復(fù)合材料機(jī)械連接結(jié)構(gòu)。

(3) 為了提高直接約束法的計(jì)算精度，可通過細(xì)分接觸面網(wǎng)格的方式獲取更為準(zhǔn)確的位移約束。

(4) 采用罰函數(shù)方法時(shí)需合理選擇罰函數(shù)使其滿足接觸剛度與接觸體剛度相匹配，以避免計(jì)算不收斂。

(5) 對(duì)于復(fù)合材料機(jī)械連接結(jié)構(gòu)中的接觸問題，推薦采用有限滑移接觸檢測(cè)技術(shù)和罰函數(shù)法接觸約束施加方法，兼顧建模效率及計(jì)算精度。

[1] Chang F K, Scott R A, Springer G S. Failure strength of nonlinearly elastic composite laminates containing a pin loaded hole[J]. Journal of Composite Materials, 1984, 18(5): 464-477.

[2] Hahn H T, Tsai S W. Nonlinear elastic behavior of unidirectional composite laminae[J]. Journal of Composite Materials, 1973, 7(1): 102-118.

[3] 陳浩然, 息志臣. 復(fù)合材料機(jī)械連接件非線性接觸應(yīng)力分析[J]. 航空學(xué)報(bào), 1990, 11(7): 410-414.

Chen Haoran, Xi zhichen. Nonlinear contact stress analysis of composite bolted joints[J]. Acta Aeronautica et Astronautica Sinica, 1990, 11(7): 410-414.(in Chinese)

[4] 張寶生, 陳家慶, 蔣力培, 等. Marc在接觸分析中的應(yīng)用[J]. 北京石油化工學(xué)院學(xué)報(bào), 2003, 11(2): 38-41.

Zhang Baosheng, Chen Jiaqing, Jiang Lipei, et al. Application of advanced nonlinear finite element analysis software Marc in contact analysis[J]. Journal of Beijing Institute of Petro-chemical Technology, 2003, 11(2): 38-41.(in Chinese)

[5] 龔思楚, 張憲政, 梅李霞, 等. 基于ABAQUS接觸算法結(jié)構(gòu)強(qiáng)度分析[J]. 教練機(jī), 2016, (4): 28-31.

Gong Sichu, Zhang Xianzheng, Mei Lixia, et al. Structure strength analysis based on ABAQUS contact arithmetics[J]. Trainer, 2016, (4): 28-31.(in Chinese)

[6] 陳向明, 柴亞南, 沈真. 自定義粘接元接觸技術(shù)在復(fù)合材料連接強(qiáng)度分析中的應(yīng)用[J]. 復(fù)合材料學(xué)報(bào), 2011, 28(6): 230-236.

Chen Xiangming, Chai Yanan, Shen Zhen. Application of contact technology based on the user-defined cohesive element in the strength analysis of composite fastened joint[J]. Acta Materiae Compositae Sinica, 2011, 28(6): 230-236.(in Chinese)

[7] 鄭潔, 李龍彬, 錢超. 復(fù)合材料機(jī)械連接層合板孔邊應(yīng)力集中研究[J]. 航空工程進(jìn)展, 2012, 3(4), 438-441.

Zheng Jie, Li Longbin, Qian Chao. Study on hole edge stress concentration of composite laminates[J]. Advances in Aeronautical Science and Engineering, 2012, 3(4), 438-441.(in Chinese)

[8] Zhou Y H, Nezhad H Y, Hou C, et al. A three dimensional implicit finite element damage model and its application to single-lap multi-bolt composite joints with variable clearance[J]. Composite Structures, 2015, 131, 1060-1072.

[9] Zhou Y H, Nezhad H Y, McCarthy M A, et al. A study of intra-laminar damage in double-lap, multi-bolt, composite joints with variable clearance using continuum damage mechanics[J]. Composite Structures, 2014, 116, 441-452.

[10] Sang H Lee. MSC/NASTRAN handbook for nonlinear analysis: based on Version 67[M]. USA: The Macneal-Schwendler Corporation, 1992.

[11] Dassault Systemes. ABAQUS analysis user’s manual[M]. USA: Dassault Systemes Simulia Corporation, 2012.

[12] MSC. Software Corporation. MARC 2007R1 user’s manual, Volume A(Theory and user information)[M]. USA: MSC Software Corporation, 2007.

[13] Bath K J, Bouzinov P A. On the constraint function method for contact problems[J]. Computer and Structures, 1997, 64(5/6): 1069-1085.

[14] Pantano A, Averill R C. A penalty-based finite element interface technology[J]. Computers and Structures, 2002, 80(22): 1725-1748.

[15] Zavarise G, Lorenzis L D. A modified node-to-segment algorithm passing the contact patch test[J]. International Journal for Numerical Methods in Engineering, 2009, 79(4): 379-416.

[16] Mayer M H, Gaul L. Segment-to-segment contact elements for modelling joint interfaces in finite element analysis[J]. Mechanical Systems and Signal Processing, 2007, 21(2): 724-734.

[17] Solberg J M, Papadopoulos P. An analysis of dual formulations for the finite element solution of two-body contact problems[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(25/26): 2734-2780.

[18] Wang W C, Sheu Y M. Stress analysis of bolted joints in CFRP laminates by half-fringe birefringent-coating technique[J]. Composite Structures, 1996, 34(1): 91-100.

[19] Hyer M W, Liu D. Stresses in a quasi-isotropic pin-loaded connector using photoelasticity[J]. Experimental Mechanics, 1984, 24(1): 48-53.