趙金娟　王世軍　楊　超　王詩義　楊慧新

西安理工大學(xué)，西安，710048

?

基于橫觀各向同性假定的固定結(jié)合部本構(gòu)關(guān)系及有限元模型

趙金娟王世軍楊超王詩義楊慧新

西安理工大學(xué)，西安，710048

提出了一種固定結(jié)合部的有限元建模方法。該方法將接觸層等效成均質(zhì)、虛擬的橫觀各向同性材料。根據(jù)結(jié)合面的變形和載荷關(guān)系，推導(dǎo)了虛擬材料本構(gòu)矩陣中5個獨(dú)立變量與法向應(yīng)力之間的非線性關(guān)系。將這個本構(gòu)關(guān)系引入到有限元分析中用以建立機(jī)器的固定聯(lián)接面模型。分析結(jié)果與試驗結(jié)果的對比證實(shí)了該方法的有效性。

固定連接；虛擬材料；有限元；本構(gòu)關(guān)系

0　引言

結(jié)合面對結(jié)構(gòu)的靜態(tài)和動態(tài)特性有顯著的影響。通常認(rèn)為機(jī)床中有30%～50%的剛度、90%的阻尼來自于結(jié)合面[1]。機(jī)械結(jié)合面的建模是機(jī)械結(jié)構(gòu)整機(jī)性能分析的關(guān)鍵[2-5]。

在基于有限元法的機(jī)械結(jié)構(gòu)整機(jī)性能分析中，結(jié)合面主要采用彈簧-阻尼單元建模[6]。使用這種單元需要根據(jù)單元附屬的接觸面積求得單元的剛度和阻尼。如果接觸表面的實(shí)體單元網(wǎng)格不規(guī)則，則每一個彈簧-阻尼單元的附屬面積都需要單獨(dú)計算，工作量很大。在有限元軟件中，這種節(jié)點(diǎn)對節(jié)點(diǎn)的彈簧-阻尼單元不能像實(shí)體單元那樣自動生成，需要手工操作，逐個生成，過程比較繁瑣。此外，這種單元相互之間沒有耦合關(guān)系，不能反映結(jié)合面相鄰區(qū)域之間的影響。

地質(zhì)力學(xué)中，如果巖石的節(jié)理和結(jié)構(gòu)之間的土壤層非常厚，在有限元分析時可以使用等參的薄層實(shí)體單元建模。這種單元的優(yōu)點(diǎn)是在有限元軟件中可以自動生成，也不需要計算彈簧-阻尼單元的附屬面積，但是在機(jī)械結(jié)構(gòu)的固定連接面中，接觸層的厚度并不確定，很難利用這種單元建立機(jī)械結(jié)合面模型[7-12]。

文獻(xiàn)[13]給出了一種利用薄層單元建立結(jié)合面模型的方法，文中假定薄層單元的厚度接近于0，結(jié)合面被看作一種虛擬的正交各向異性材料并且給出了這種虛擬材料的正交各向異性的本構(gòu)矩陣：

[σxxσyyσzzσxyσyzσzx]T=

diag(0，0，E33，0，G，G) [εxxεyyεzzεxyεyzεzx]T

(1)

這里的E33是虛擬材料沿結(jié)合面法線方向的彈性模量，G是虛擬材料的剪切模量。在這個對角形式的本構(gòu)關(guān)系中，有3個對角元素等于0，這使得矩陣是奇異的。材料力學(xué)中，材料本構(gòu)矩陣的對角元素必須大于零[14]，因此式(1)的本構(gòu)關(guān)系并不能直接用于結(jié)合面的建模。此外，接近于0的厚度在結(jié)合面幾何建模時也很難操作。

本文提出了一種基于虛擬材料的具有有限厚度的固定結(jié)合面建模方法。該方法將結(jié)合面及其附近區(qū)域等效成均勻的橫觀各向同性區(qū)域并用八節(jié)點(diǎn)的六面體實(shí)體單元模擬。接觸層的實(shí)體單元可以在有限元軟件中自動生成，不需要計算單元的附屬面積，同時實(shí)體單元也能很好地反映結(jié)合面相鄰區(qū)域的相互作用。

1　虛擬材料的本構(gòu)關(guān)系和接觸單元的剛度矩陣

如果材料在平行于名義接觸平面oxy的方向上具有相同的機(jī)械性質(zhì)，則有下面的應(yīng)變-應(yīng)力關(guān)系[14]：

(2)

式中，Ex、Ey、Ez分別為x、y、z方向上的彈性模量；μij為單獨(dú)在j方向作用正應(yīng)力σj而無其他應(yīng)力分量時，i方向應(yīng)變與j方向應(yīng)變之比的負(fù)值，即泊松比，i, j=x, y, z; Gxy、Gyz為Gxz剪切模量，下標(biāo)的第一個字母表示法線方向，第二個字母表示剪切變形的方向。

式(2)是一個包含獨(dú)立參數(shù)Ex、Ez、μxy、μxz和Gxz的橫觀各向同性材料的本構(gòu)方程。

圖1a所示的包含接觸界面的矩形微體中，接觸層平行于oxy平面。假定微體的厚度dz遠(yuǎn)大于接觸層的厚度，外表面受到均布的法向和切向載荷。根據(jù)微體的應(yīng)變-應(yīng)力關(guān)系，可以求得圖1b所示的虛擬彈性體(與圖1a中微體等效、均一的橫觀各向同性的虛擬彈性體)的5個獨(dú)立參數(shù)。

(a)微體　　　　　　(b)虛擬彈性體圖1　包含接觸層的微體和等效的虛擬彈性體

1.1彈性模量Ex

圖1a中，在微體表面施加x方向的法向載荷σx，x方向產(chǎn)生應(yīng)變εx。微體的厚度比接觸層的真實(shí)厚度大很多，可以認(rèn)為圖1a中包含結(jié)合面的微體的應(yīng)變εx等于圖1b中不包含結(jié)合面的微體的應(yīng)變εx，接觸層的存在與否，并不影響σx與εx之間的關(guān)系。σx與εx仍然保持基體材料的應(yīng)力-應(yīng)變關(guān)系。圖1b中,虛擬彈性體的彈性模量Ex等于基體材料的彈性模量E。

1.2彈性模量Ez

(3)

(4)

由于無法確定準(zhǔn)確的接觸層厚度，所以接觸層的應(yīng)變ε″z也無法計算。本文按下式定義接觸層的應(yīng)變：

ε″z=λn/dz

(5)

其中，λn為接觸層在法向載荷σz下的法向變形，λn-σz的關(guān)系可以通過試驗或者理論分析獲得。根據(jù)這個定義，接觸層的應(yīng)變與微體的厚度dz相關(guān)。隨著微體厚度的變化，接觸層的應(yīng)變ε″z也隨著改變。試驗研究已經(jīng)顯示接觸層的法向變形是法向載荷的冪函數(shù)[15-16]：

(6)

其中，c和m是系數(shù)，可通過擬合試驗數(shù)據(jù)確定。因此，接觸層的應(yīng)變ε″z可用寫成下面的形式：

(7)

在施加法向載荷σz后，微體將產(chǎn)生如下的法向應(yīng)變：

(8)

由于

σz=Ezεz

(9)

所以微體在z方向的彈性模量

(10)

1.3泊松比μxy

泊松比μxy反映了x方向法向應(yīng)變εx與y方向法向應(yīng)變εy之間的耦合關(guān)系。本文認(rèn)為μxy與基體材料的泊松比μ相同：

μxy=μ

(11)

1.4泊松比μxz

(12)

(13)

(14)

接觸層的應(yīng)變ε″z與式(5)、式(7) 相同：

(15)

這樣，虛擬材料的泊松比

(16)

1.5剪切模量Gxz

Gxz是虛擬彈性體在oxz平面內(nèi)的剪切模量。圖1中的微體被假定是橫觀各向同性的，微體在oxz面內(nèi)的剪切剛度與oyz面內(nèi)的剪切剛度相等：

Gxz=Gyz

(17)

(18)

不包含接觸層時，微體在x方向的剪應(yīng)變

(19)

接觸層的剪應(yīng)變

γ″xz=λxz/dz

(20)

其中，λxz是接觸層在x方向的剪切變形。類似于接觸層的法向變形λn，λxz也能夠從結(jié)合面試驗獲得：

(21)

其中，ατ、βτ都是與法向載荷σz相關(guān)的系數(shù)，與式(7)中的c和m類似。根據(jù)式(18)～式(21)，可以得到

(22)

即

(23)

這樣，oxz面內(nèi)的剪切模量

(24)

最終，如果通過試驗獲得c、m、ατ、βτ，就可以求得虛擬材料的5個獨(dú)立變量：

(25)

1.6接觸單元的剛度矩陣

有限元分析中，八節(jié)點(diǎn)的三維實(shí)體等參元是一種經(jīng)常使用的實(shí)體單元，大多數(shù)商業(yè)軟件也提供這種單元，它的剛度矩陣為[17]

(26)

式中，B、D、J分別為應(yīng)變矩陣、彈性矩陣和雅可比行列式；wξ,i、wη,j、wζ,k為高斯積分的權(quán)系數(shù)。

由于式(2)中的系數(shù)矩陣的逆即為彈性矩陣D，根據(jù)結(jié)合面的試驗數(shù)據(jù)，就可以得到虛擬材料的彈性矩陣。將這個彈性矩陣代入式(26)，即可得到薄層單元的剛度矩陣。接觸單元的形狀不同于通常的三維實(shí)體單元，單元的2個表面必須平行于名義的接觸界面，以便保持單元厚度的一致。接觸層的真實(shí)厚度非常小，采用太厚的接觸單元會顯著改變零件中的應(yīng)力分布。在前述5個獨(dú)立變量的分析中，假定微體的厚度必須遠(yuǎn)大于真實(shí)接觸層的厚度，因此接觸單元的厚度也必須遠(yuǎn)大于真實(shí)接觸層的厚度。在微體的應(yīng)變分析中，所有的應(yīng)變都假定是小應(yīng)變。在接觸層變形相同的條件下，單元厚度越小，單元的應(yīng)變越大，這使得單元的厚度不能太小。文獻(xiàn)[12]認(rèn)為等參單元的邊長比超過1000也不會產(chǎn)生顯著的數(shù)值誤差，因而這種等參元用作薄層單元模擬接觸層的性質(zhì)是可行的。

2　本構(gòu)關(guān)系和接觸單元的試驗驗證

文獻(xiàn)[18]通過試驗研究了結(jié)合面切向變形與切向載荷的關(guān)系。試驗結(jié)果顯示，結(jié)合面出現(xiàn)宏觀滑動之前，經(jīng)歷了彈性變形、塑性變形和微滑移三個階段。結(jié)合面只在彈性變形階段才表現(xiàn)出線性的切向剛度。在彈塑性過渡階段和微滑移階段，結(jié)合面表現(xiàn)出非線性的切向剛度。圖2為文獻(xiàn)[18]中試驗裝置的示意圖。2個100 mm×50 mm×18 mm的鋼制磨削試樣之間的接觸面為研究的結(jié)合面。結(jié)合面的法向載荷是試樣A的自重以及上面的配重。試樣B固定在基礎(chǔ)梁上，基礎(chǔ)梁通過減振器固定在工作臺上，這樣能夠隔絕來自地面的振動。結(jié)合面的切向載荷施加在試樣A的兩側(cè)并且載荷作用點(diǎn)靠近結(jié)合面以避免在結(jié)合面上產(chǎn)生力矩。

圖2　切向變形的試驗裝置

結(jié)合面的切向變形通過2個分辨率為12.5 nm的電容傳感器測量。如圖2所示，傳感器SB放在試樣B的一端，測量試樣B的切向位移，傳感器SA用來測量試樣A的切向位移。如果接觸層的法向應(yīng)力和切向應(yīng)力是均勻的，那么2個傳感器測量結(jié)果的差值就是接觸層的切向變形。

試驗裝置的有限元模型包含1008個單元，其中，144個單元是本文提出的接觸單元。這個模型中，接觸層的厚度是100 μm，由一層矩形六面體接觸單元劃分而成。與圖2相同，試樣B的下表面被固定，20 kPa的壓力施加在試樣A的上表面，每側(cè)的切向力是3 N。試樣的材料參數(shù)如下：彈性模量E=200 GPa，泊松比μ=0.3，密度ρ=7850 kg/m3。接觸層的試驗數(shù)據(jù)采用文獻(xiàn)[16]中的數(shù)據(jù)：

(27)

其中，pn(pn≤2.5 MPa)為結(jié)合面的法向壓力，相當(dāng)于式(25)中的σz；knj、ktj分別為結(jié)合面單位面積上的法向剛度和切向剛度，對照式(6)和式(21)可知

(28)

式(28)與式(27)在形式上相同，對照兩式即可確定參數(shù)c、m、ατ、βτ，進(jìn)而可以確定式(25)中虛擬材料的5個獨(dú)立變量。

式(25)中，有些材料參數(shù)與接觸層的法向載荷之間是非線性關(guān)系，這使得有限元模型是非線性的。如果接觸單元的法向應(yīng)力不一樣，則接觸單元的材料參數(shù)也不一樣。計算時，首先施加法向壓力pn，當(dāng)法向壓力pn增加到20 kPa時，開始施加切向載荷pt。計算獲得的剪切變形與剪切載荷之間的關(guān)系顯示在圖3中。作為比較，文獻(xiàn)[18]的試驗結(jié)果也顯示在圖3中。對比兩條曲線可以發(fā)現(xiàn)，計算結(jié)果與試驗結(jié)果是一致的。文獻(xiàn)[18]的試驗中，法向載荷只有20 kPa，當(dāng)切向載荷pt超過4.7 kPa時，切向變形開始進(jìn)入塑性變形階段并出現(xiàn)明顯的宏觀滑移；切向載荷小于4.7 kPa時，宏觀滑移較小，切向變形與切向載荷大體呈線性關(guān)系，但是比計算結(jié)果略大。式(27)適用于最大法向載荷為2.5 MPa的情況，遠(yuǎn)大于文獻(xiàn)[18]中20 kPa的法向載荷，式(27)的適用條件與文獻(xiàn)[18]中的試驗條件并不完全相同。由于結(jié)合面的非線性性質(zhì)[1]，在切向載荷相同的條件下，法向載荷大的結(jié)合面，其切向剛度也大。因此圖3中切向變形的試驗結(jié)果比計算結(jié)果略大是合理的。

圖3　切向變形的試驗結(jié)果和有限元結(jié)果比較

當(dāng)法向壓力pn=20 kPa，切向載荷pt=4.7 kPa時，接觸層沿x方向的剪切應(yīng)力如圖4所示。圖4中，接觸層的剪切應(yīng)力分布不均勻，應(yīng)力變化達(dá)到9.5%，這使得接觸層的剪切變形也不一致。出現(xiàn)這種現(xiàn)象的原因在于切向載荷不是均勻施加在接觸層上的。這種不一致的剪切變形，使得傳感器A和B的測量結(jié)果不能夠準(zhǔn)確地反映接觸層的切向變形，這也是結(jié)合部切向特性的試驗研究中的主要困難之一。

圖4　接觸層x方向的剪應(yīng)力分布

3　結(jié)語

本文提出的利用等參實(shí)體單元建立固定結(jié)合面模型的方法中，結(jié)合面被等效為一個虛擬的橫觀各向同性實(shí)體，實(shí)體單元的材料特性通過結(jié)合面的試驗結(jié)果獲得。常用的有限元軟件都有本文方法所使用的材料本構(gòu)關(guān)系和單元類型，因此本文的方法能夠方便地應(yīng)用在多數(shù)商業(yè)性的有限元軟件中。實(shí)體單元可以通過軟件快速自動生成，從而給結(jié)合面的有限元建模帶來便利。

本文采用試驗獲得的結(jié)合面變形-載荷關(guān)系建立了虛擬材料的本構(gòu)矩陣。事實(shí)上，本構(gòu)關(guān)系的建立也可以利用其他方法獲得的結(jié)合面的變形-載荷關(guān)系，例如基于粗糙表面輪廓的有限元分析，分形接觸理論以及統(tǒng)計接觸理論等，本文提出的結(jié)合面建模方法具有通用性。

[1]趙宏林, 丁慶新, 曾鳴，等. 機(jī)床結(jié)合部特性的理論解析及應(yīng)用[J]. 機(jī)械工程學(xué)報, 2008, 44(12):208-213.

ZhaoHonglin,DingQingxin,ZengMing,etal.TheoreticAnalysisonandApplicationofBehaviorsofMachineToolJoints[J].JournalofMechanicalEngineering, 2008, 44(12): 208-213.

[2]石坤, 宋俐, 師俊平. 機(jī)械結(jié)合部等效材料參數(shù)建立與試驗[J]. 農(nóng)業(yè)機(jī)械學(xué)報, 2014, 44(6): 297-301.

ShiKun,SongLi,ShiJunping.EstablishmentandExperimentofMechanicalJointEquivalentMaterialProperties[J].TransactionsoftheChineseSocietyforAgriculturalMachinery, 2014, 44(6): 297-301.

[3]黎定仕, 張以都, 王鵬. 基于結(jié)構(gòu)阻尼的機(jī)械結(jié)合部動力學(xué)模型研究[J]. 振動與沖擊, 2010, 29(8): 204-208.

LiDingshi,ZhangYidu，WangPeng.DynamicModelofMachineJointsBasedonStructuralDamping[J].JournalofVibrationandShock, 2010, 29(8): 204-208.

[4]米良, 殷國富, 孫明楠, 等. 結(jié)合部動力學(xué)特性的立柱-主軸系統(tǒng)動力學(xué)模型研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報, 2011, 42(12): 202-207.

MiLiang,YinGuofu,SunMingnan,etal.ColumnspindleSystemDynamicModelBasedonDynamicCharacteristicsofJoints[J].TransactionsoftheChineseSocietyforAgriculturalMachinery, 2011, 42(12):202-207.

[5]劉恒, 劉意, 王為民. 接觸界面法向剛度等效的新方法[J]. 機(jī)械工程學(xué)報, 2011, 47(17): 37-43.LiuHeng,LiuYi,WangWeimin.NewEquivalentMethodforNormalStiffnessofContactInterface[J].JournalofMechanicalEngineering, 2011, 47(17): 37-43.

[6]TianHongliang,LiBin,LiuHongqi,etal.ANewMethodofVirtualMaterialHypothesis-basedDynamicModelingonFixedJointInterfaceinMachineTools[J].InternationalJournalofMachineTools&Manufacture, 2011, 51: 239-249.

[7]AhmadianH,EbrahimiM,MottersheadJ,etal.IdentificationofBoltedJointsInterfaceModels[C]//ProceedingsofISMA2002:NoiseandVibrationEngineering.Leuven,Belgium, 2002: 1741-1747.

[8]AhmadianH,JalaliH,MottersheadJ,etal.DynamicModelingofSpotWeldsUsingThinLayerInterfaceTheory[C]//ProceedingsoftheTenthInt.CongressonSoundandVibrationICSV10.Stockholm,Sweden, 2003: 3439-3446.

[9]AhmadianH,MottersheadJ,JamesS,etal.ModellingandUpdatingofLargeSurface-to-surfaceJointsintheAWE-MACEStructure[J].MechanicalSystemsandSignalProcessing, 2006, 20: 868-880.

[10]DesaiC,ZamanM,LightnerJ,etal.Thin-layerElementforInterfacesandJoints[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 1984, 8: 19-43.

[11]SharmaK,DesaiC.AnalysisandImplementationofThin-layerElementforInterfacesandJoints[J].JournalofEngineeringMechanics， 1992, 118: 2442-2462.

[12]PandeG,SharmaK.OnJoint/interfaceElementsandAssociatedProblemsofIll-conditioning[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 1979, 3: 293-300.

[13]BogradS,ReussP,ASchmidt,etal.ModelingtheDynamicsofMechanicalJoints[J].MechanicalSystemsandSignalProcessing, 2011, 25: 2801-2826.

[14]劉新東, 劉偉. 復(fù)合材料力學(xué)基礎(chǔ)[M]. 西安: 西北工業(yè)大學(xué)出版社, 2010.

[15]FuWP,HuangYM,ZhangXL,etal.ExperimentalInvestigationofDynamicNormalCharacteristicsofMachinedJointSurfaces[J].JournalofVibrationandAcoustics,TransactionsoftheASME, 2000, 122: 393-398.

[16]石坤, 宋俐, 師俊平, 等. 一種新的機(jī)械結(jié)構(gòu)結(jié)合部特性分析方法[J]. 機(jī)械工程學(xué)報, 2013, 49(1): 142-147.

ShiKun,SongLi,ShiJunping,etal.ANewMethodforCharacteristicAnalysisoftheMechanicalStructureJoint[J].JournalofMechanicalEngineering, 2013, 49(1): 142-147.

[17]秦太驗, 徐春暉, 周喆. 有限元法及其應(yīng)用[M]. 北京: 中國農(nóng)業(yè)大學(xué)出版社, 2011.

[18]NiJun,ZhuZhenqi.ExperimentalStudyofTangentialMicroDeflectionofInterfaceofMachinedSurfaces[J].JournalofManufacturingScienceandEngineering, 2001, 123: 365-367.

(編輯張洋)

A Constitutive Law Based on Transverse Isotropic Hypothesis and Finite Element Model of Fixed Joints

Zhao JinjuanWang ShijunYang ChaoWang ShiyiYang Huixin

Xi’an University of Technology，Xi’an，710048

A finite element method modeling fixed joints was presented. In the method, the contact layers of the fixed joints were regarded as a homogenous virtual transverse isotropic material. According to the relations among deformations and loads of joints, the nonlinear relations among 5 independent variables and normal stress in constitutive matrix were deduced, and then the relations were introduced into finite element analysis to model fixed joints. The analytical results and experimental results of the joints were compared and the validity of the method was confirmed.

fixed joint； virtual material； finite element； constitutive law

趙金娟，女，1974年生。西安理工大學(xué)印刷包裝工程學(xué)院講師。主要研究方向為機(jī)械接觸特性及建模方法。發(fā)表論文15篇。王世軍(通信作者)，男，1967年生。西安理工大學(xué)機(jī)械與精密儀器學(xué)院副教授。楊超，男，1988年生。西安理工大學(xué)機(jī)械與精密儀器學(xué)院碩士研究生。王詩義，男，1989年生。西安理工大學(xué)機(jī)械與精密儀器學(xué)院碩士研究生。楊慧新，男，1992年生。西安理工大學(xué)機(jī)械與精密儀器學(xué)院碩士研究生。

2015-05-29

國家重點(diǎn)基礎(chǔ)研究計劃(973計劃)資助項目(2009CB724406)；陜西省科技廳科技統(tǒng)籌創(chuàng)新工程重點(diǎn)實(shí)驗室資助項目(2014SZS10-P05)

TH123

10.3969/j.issn.1004-132X.2016.08.003