基于橫觀各向同性假定的固定結合部本構關系及有限元模型

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

西安理工大學，西安，710048

?

基于橫觀各向同性假定的固定結合部本構關系及有限元模型

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

西安理工大學，西安，710048

提出了一種固定結合部的有限元建模方法。該方法將接觸層等效成均質、虛擬的橫觀各向同性材料。根據結合面的變形和載荷關系，推導了虛擬材料本構矩陣中5個獨立變量與法向應力之間的非線性關系。將這個本構關系引入到有限元分析中用以建立機器的固定聯接面模型。分析結果與試驗結果的對比證實了該方法的有效性。

固定連接；虛擬材料；有限元；本構關系

0　引言

結合面對結構的靜態和動態特性有顯著的影響。通常認為機床中有30%～50%的剛度、90%的阻尼來自于結合面[1]。機械結合面的建模是機械結構整機性能分析的關鍵[2-5]。

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

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

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

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

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

(1)

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

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

1　虛擬材料的本構關系和接觸單元的剛度矩陣

如果材料在平行于名義接觸平面oxy的方向上具有相同的機械性質，則有下面的應變-應力關系[14]：

(2)

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

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

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

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

1.1彈性模量Ex

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

1.2彈性模量Ez

(3)

(4)

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

ε″z=λn/dz

(5)

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

(6)

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

(7)

在施加法向載荷σz后，微體將產生如下的法向應變：

(8)

由于

σz=Ezεz

(9)

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

(10)

1.3泊松比μxy

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

μxy=μ

(11)

1.4泊松比μxz

(12)

(13)

(14)

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

(15)

這樣，虛擬材料的泊松比

(16)

1.5剪切模量Gxz

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

Gxz=Gyz

(17)

(18)

不包含接觸層時，微體在x方向的剪應變

(19)

接觸層的剪應變

γ″xz=λxz/dz

(20)

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

(21)

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

(22)

即

(23)

這樣，oxz面內的剪切模量

(24)

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

(25)

1.6接觸單元的剛度矩陣

有限元分析中，八節點的三維實體等參元是一種經常使用的實體單元，大多數商業軟件也提供這種單元，它的剛度矩陣為[17]

(26)

式中，B、D、J分別為應變矩陣、彈性矩陣和雅可比行列式；wξ,i、wη,j、wζ,k為高斯積分的權系數。

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

2　本構關系和接觸單元的試驗驗證

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

圖2　切向變形的試驗裝置

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

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

(27)

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

(28)

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

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

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

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

圖4　接觸層x方向的剪應力分布

3　結語

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

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

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

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

[2]石坤, 宋俐, 師俊平. 機械結合部等效材料參數建立與試驗[J]. 農業機械學報, 2014, 44(6): 297-301.

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

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

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

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

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

[5]劉恒, 劉意, 王為民. 接觸界面法向剛度等效的新方法[J]. 機械工程學報, 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]劉新東, 劉偉. 復合材料力學基礎[M]. 西安: 西北工業大學出版社, 2010.

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

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

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

[17]秦太驗, 徐春暉, 周喆. 有限元法及其應用[M]. 北京: 中國農業大學出版社, 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年生。西安理工大學印刷包裝工程學院講師。主要研究方向為機械接觸特性及建模方法。發表論文15篇。王世軍(通信作者)，男，1967年生。西安理工大學機械與精密儀器學院副教授。楊超，男，1988年生。西安理工大學機械與精密儀器學院碩士研究生。王詩義，男，1989年生。西安理工大學機械與精密儀器學院碩士研究生。楊慧新，男，1992年生。西安理工大學機械與精密儀器學院碩士研究生。

2015-05-29

國家重點基礎研究計劃(973計劃)資助項目(2009CB724406)；陜西省科技廳科技統籌創新工程重點實驗室資助項目(2014SZS10-P05)

TH123

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