陳曉霞，程 凱

（1.天津工業(yè)大學機械工程學院，天津 300387；2.天津工業(yè)大學天津市現代機電裝備技術重點實驗室，天津300387）

基于三維實體有限元模型的諧波齒輪裝配狀態(tài)應力分析

陳曉霞1，2，程 凱1

（1.天津工業(yè)大學機械工程學院，天津 300387；2.天津工業(yè)大學天津市現代機電裝備技術重點實驗室，天津300387）

為準確地反映諧波柔輪的齒廓及齒根的應力狀態(tài)，提出了基于真實齒廓信息的柔輪三維實體的建模方法.以漸開線齒形為例，利用APDL語言建立包括漸開線齒廓、齒厚和齒根倒圓等參數的實體模型.定義波發(fā)生器與柔輪筒內表面間接觸，求解獲得波發(fā)生器作用下的柔輪應力，與基于圓環(huán)理論的理論值比較分析；為更真實反映柔輪在裝配狀態(tài)下的應力分布，構建柔輪三維實體模型，求解柔輪應力.通過路徑定義，提取與理論計算截面相對應的中截面上的應力，與平面齒圈的有限元結果和理論值比較分析.研究發(fā)現：平面齒圈有限元實體模型中，沿周向應力波動均值與理論值基本吻合，驗證了文中的有限元實體建模方法的有效性.三維實體模型周向應力比平面齒圈周向應力大，且最大周向應力發(fā)生在頂層.

諧波齒輪；柔輪；有限元分析；三維實體模型；裝配狀態(tài)

諧波齒輪傳動是基于柔輪的彈性變形來傳遞運動的.諧波齒輪的工作壽命主要取決于柔輪的工作壽命，因此柔輪的強度研究是諧波齒輪傳動中的重要課題，而確定柔輪中的應力大小和分布規(guī)律是強度研究的關鍵［1］.

目前，許多研究者基于實體有限元模型來研究柔輪的應力和變形［2-3］.在使用ANSYS過程中，為了計算方便，研究者們都采用了不同的簡化方法［4-5］.王知行等［6］把柔輪的輪齒均勻地抹平，看作光滑變厚度的薄壁殼；韓建華等［7］把齒圈視作當量厚度的光殼，分析了二點式波發(fā)生器作用下柔輪的應力分布規(guī)律，得出了柔輪結構參數對應力分布的影響和柔輪發(fā)生疲勞破壞的可能截面；Chen等［8-9］用殼單元和梁單元相結合的柔輪有限元模型計算雙圓盤波發(fā)生器作用下柔輪變形和應力，得到齒圈部位的中性層、齒根、齒尖處的變形；付軍峰等［10］建立了諧波齒輪傳動中杯形柔輪與波發(fā)生器接觸分析的有限元模型，得到了柔輪殼體空載時的應力分布狀況.

劉文芝等［11］在嚙合齒對接觸邊界的節(jié)點之間建立面-面接觸單元，用三維彈性接觸問題有限元法較全面地計算和分析了承載柔輪齒圈和筒體的應力和分布規(guī)律、輪齒變形及影響機理；吳偉國等［12］在ANSYS中建立帶有整圈輪齒的剛輪和柔輪模型，利用其接觸分析要素在嚙合齒面間定義面-面接觸對，使得各對輪齒在嚙合時產生作用力，從而提出了一種與實際較接近的諧波傳動有限元模型建模方法；董慧敏等［13-15］將柔輪輪齒簡化為等厚度的殼體，用經驗公式和實驗結果計算得到了柔輪在不同負載條件下齒圈中面不同截面上的變形情況.

本文針對現有諧波傳動有限元分析模型不夠精細這一問題，利用APDL語言建立包括柔輪齒厚、齒根倒圓半徑和齒槽壁厚等參數的參數化實體模型，并結合柔輪的漸開線方程，參數化設計了柔輪的漸開線齒廓，生成平面齒圈，添加平面波發(fā)生器.在平面齒圈徑向上定義頂層、中層、底層3個路徑.定義波發(fā)生器與柔輪筒內表面之間為線-線接觸，求解獲得波發(fā)生器作用下的柔輪應力，與基于圓環(huán)理論的理論值和設計規(guī)范值比較分析.將平面齒圈和波發(fā)生器軸向拉伸，再添加柔輪筒，生成柔輪的三維模型.在柔輪三維實體模型徑向上定義頂層、中層、底層3個路徑，軸向上定義前截面、中截面、后截面3個路徑.定義波發(fā)生器與柔輪筒內表面之間為面-面接觸，求解獲得波發(fā)生器作用下的柔輪應力.提取三維模型中截面上不同層面上的應力與平面齒圈相應層上的有限元值、理論值和設計規(guī)范值進行比較分析.

1 裝配狀態(tài)下的應力計算

1.1 基于等效圓環(huán)理論的應力計算

本次研究采用四滾輪波發(fā)生器，以杯形諧波齒輪為例，剛輪和柔輪均采用變位的漸開線齒形.

按照圓環(huán)理論，考慮柔輪筒對齒圈的加強作用，建立具有等效剛度的等效圓環(huán)模型.利用等效圓環(huán)模型可以得到空載傳動狀態(tài)的圓環(huán)變形和應力計算公式，為嚙合分析和應力強度計算提供依據.通過理論計算可以得到柔輪任意截面上的周向應力［16］：

式中：φ為波發(fā)生器長軸相對柔輪未變形端的轉角；β為滾輪安裝軸的偏置角度；w0為最大徑向變形量；E為材料的彈性模量；h為圓環(huán)的厚度，h=1.671/3s［10］，s為柔輪齒圈處的壁厚；rm為變形前柔輪中性圓曲線的半徑.

1.2 設計規(guī)范中的應力計算

在設計規(guī)范中，為了簡化強度計算，通常把柔輪簡化為一個光滑圓柱殼體進行應力分析，然后再根據試驗結果進行適當的修正.可求得，柔輪的周向應力最大值為［17］

式中：Krt為柔輪彎曲剛度引起的應力增大系數，設s*= s/m，Krt=（1+s*）/s*；KM為受載柔輪畸變引起的應力增大系數；Kd為動載荷系數；Cσ為正應力系數.

2 柔輪三維實體模型的建立

2.1 主要幾何參數

為了實現模型的參數化設計，本文模型的建立、分析、后處理等過程采用APDL語言.建模時，為計算方便，空載傳動狀態(tài)下的柔輪齒圈采用1/4柔輪模型對稱約束；在β=25°處建立四滾輪波發(fā)生器并添加接觸定義，波發(fā)生器作用在齒圈的中部；柔輪齒根過渡圓角處網格劃分較密.

柔輪輪齒示意圖如圖1所示，選用有限元軟件ANSYS進行建模和分析.

選取柔輪輪齒的結構參數如表1所示.

2.2 平面齒圈的建立

當齒輪用滾刀或插刀切削時，漸開線齒廓表達式為［9］

圖1 柔輪輪齒示意圖Fig.1 Diagram of flexpline tooth profile

表1 柔輪輪齒的參數Tab.1 Parameters of flexible tooth

式中：r1為柔輪分度圓半徑；u為系刀具沿被切齒輪的分度圓作無滑動的滾動時之轉角；α0為基準齒形角；θ為齒輪分度圓齒厚所對中心角的1/2.

柔輪采用滾動加工時，其漸開線齒廓的幾何尺寸計算公式如表2所示.

表2 漸開線柔輪齒廓的幾何尺寸計算公式Tab.2 Calculation formulas for geometric dimensions of involute tooth profile of flexpline

根據圖1的柔輪輪齒簡圖和表1、表2中的柔輪參數，利用APDL語言，確定柔輪的齒頂圓和齒根圓的位置.再結合柔輪的漸開線方程式（3）確定漸開線齒廓的點的坐標.采用點-線-面的方法，生成半個齒，再將半個齒鏡像生成一個完整的齒，通過陣列生成1/4輪齒，如圖2所示.

圖2 平面齒圈模型Fig.2 Model of planar tooth ring

2.3 三維實體模型的建立

將平面齒圈沿軸向拉伸，生成三維實體模型.基于圖3結構模型和表3結構參數添加柔輪筒體.

圖3 柔輪結構簡圖Fig.3 Structure diagram of flexpline

表3 柔輪筒的結構參數Tab.3 Structural parameters of flexible cylinder

利于APDL語言，在網格剖分時，將齒的底面與柔輪筒上與實體齒圈接觸的部分周向均分和軸向均分一致，使兩者剖分后的節(jié)點位置一一對應.通過APDL語言中“Nummrg”命令來實現節(jié)點合并.獲得如圖4所示三維實體模型.

圖4 柔輪三維實體模型Fig.4 3D solid model of flexpline

3 實體有限元模型的應力分析

3.1 平面齒圈的應力分析

計算理論值和規(guī)范值時，選取中性圓半徑rm= 50.375 mm，選取徑向變形系數為1，則平面齒圈最大徑向位移w0=m=0.5 mm.取E=210 GPa，KM=1，Kd=1，Cσ=1.5［17］.

平面齒圈中，沿柔輪筒的徑向壁厚分為3層.齒根圓半徑處設為頂層，柔輪筒中面設為中層，柔輪筒的內壁設為底層，如圖5所示.

圖5 徑向分層示意圖Fig.5 Radial layer diagram of flexpline

定義波發(fā)生器與柔輪筒內表面之間為線-線接觸，利用ANSYS軟件對平面齒圈進行有限元計算分析，得到柔輪的應力云圖，如圖6所示.

圖6 平面齒圈應力云圖Fig.6 Stress nephogram of planar tooth ring

根據在圖5所示定義的3條路徑，提取頂層、中層、底層應力.由于中層應力較小，在-6.5～7.5 MPa之間，對柔輪影響不大，不作討論.提取平面齒圈的頂層應力與式（1）的理論值比較，如圖7所示.齒圈發(fā)生彎曲變形時，頂層受到拉應力，底層受到壓應力，提取平面齒圈的底層應力，與由式（1）計算得到的理論值的相反數比較，可以得到比較曲線如圖8所示.

圖7 平面齒圈頂層周向應力與理論值比較Fig.7 Circumferential stress distribution in top layer of planar tooth ring and theoretical result

圖8 平面齒圈底層周向應力與理論值比較Fig.8 Circumferential stress distribution in bottom layer of planar tooth ring and theoretical result

基于式（1）等效圓環(huán)理論的最大值為52.2 MPa，發(fā)生在柔輪短軸處；在柔輪與波發(fā)生器接觸處，理論值為49 MPa.由式（2）設計規(guī)范獲得最大周向彎曲應力為74.5 MPa.

通過圖7和圖8平面齒圈的應力可以發(fā)現，各層沿周向應力分布是不均勻的.應力波動幅度在柔輪的短軸、長軸和波發(fā)生器接觸處較大，在φ=45°部位最小.沿著周向應力波動的均值的變化趨勢與理論計算值基本一致.

如圖7所示，平面齒圈頂層上最大周向應力為106.2 MPa，發(fā)生在柔輪短軸處，比理論值大103.4%，比設計規(guī)范值大42.6%；在柔輪與波發(fā)生器接觸處，平面齒圈頂層周向應力為100.2 MPa，比理論值大51.1%，比設計規(guī)范值大34.5%.

如圖8所示，平面齒圈的底層上最大周向應力為66 MPa，發(fā)生在柔輪短軸處，比理論值大25.4%，比設計規(guī)范值小11.4%；在柔輪與波發(fā)生器接觸處，平面齒圈底層周向應力為61 MPa，比理論值大24.5%，比設計規(guī)范值小18.1%.

3.2 三維實體模型的應力分析

把平面齒圈軸向拉伸，生成三維實體齒圈.沿齒圈的軸向方向，把柔輪的齒寬分為圖9所示的前、中、后3個截面.

圖9 軸向分層示意圖Fig.9 Axial layer diagram of flexpline

定義波發(fā)生器與柔輪筒內表面之間為面-面接觸，利用ANSYS軟件對三維實體模型進行計算分析，求解獲得柔輪三維實體模型的應力云圖，如圖10所示.

圖10 柔輪三維實體模型應力云圖Fig.10 Stress nephogram of 3D flexpline model

平面齒圈模型與三維實體齒圈的中截面是對應的，故提取三維實體模型中截面上的應力，與平面齒圈模型計算結果進行比對分析.提取三維實體模型中截面上頂層和底層的周向應力如圖11、圖12所示.

圖11 三維實體模型中截面頂層周向應力與理論值比較Fig.11 Circumferentialstressdistributionintoplayerinmiddle section of 3D solid model and theoretical result

圖12 三維實體模型中截面底層周向應力與理論值比較Fig.12 Circumferential stress distribution in bottom layer in middle section of 3D solid model and theoretical result

如圖11所示，三維實體模型中截面頂層最大周向為170.3 MPa，發(fā)生在柔輪筒短軸處，比平面齒圈大60%，比理論值大226%，比設計規(guī)范值大129%.在波發(fā)生器與柔輪接觸處，中截面周向應力為120 MPa，比平面齒圈大20%，比理論值大145%，比設計規(guī)范值大61%.

如圖12所示，三維實體模型中截面底層最大周向為88 MPa，發(fā)生在柔輪短軸處，比平面齒圈大33%，比理論值大69%，比設計規(guī)范值大18%.在波發(fā)生器與柔輪接觸處，中截面底層周向應力為60 MPa，與平面齒圈的相等，比理論值大22%，比設計規(guī)范值小19%.

通過圖11和圖12三維實體模型的應力曲線發(fā)現，三維實體模型內沿周向應力波動的均值的變化規(guī)律與理論值基本一致.但無論在頂層，還是在底層，三維實體模型的應力計算結果都大于平面齒圈模型計算結果.這是由于平面齒圈沒有軸向應力，在三維實體模型中的應力為雙向應力，中截面在彎曲狀態(tài)下產生軸向應力，使周向應力增大.

4 結論

使用參數化建模方法，在ANSYS中利用APDL語言建立包含漸開線齒廓的柔輪有限元實體模型，通過設置波發(fā)生器與柔輪的接觸關系，分析柔輪在空載狀態(tài)下的應力，獲得更為接近柔輪實際情況的應力結果.

（1）平面齒圈有限元實體模型中，沿周向應力波動均值與理論計算值基本吻合.

（2）平面齒圈有限元實體模型的周向應力計算表明：柔輪齒圈上的周向應力在底層較小，在頂層最大；頂層齒槽底部的周向應力均比規(guī)范值大.

（3）三維實體模型中的中截面上最大周向應力也發(fā)生在頂層，但相比于平面齒圈，周向應力極值都偏大.三維實體模型中截面頂層最大周向比理論值大2倍以上，比設計規(guī)范值大1倍以上.

［1］MUSSER C W.Strain wave gearing：US，2906143［P］.1959-09-01.

［2］辛洪兵.柔輪齒圈應力的有限元分析［J］.機械科學與技術，2003，22（4）：558-559. XIN H B.Finite element analysis of the stress of flexible tooth ring［J］.Mechanical Science and Technology，2003，22（4）：558-559（in Chinese）.

［3］向國齊．諧波齒輪傳動柔輪有限元分析研究 ［D］．成都：四川大學，2005．XIANG G Q.Finite element analysis and study of flexspline in harmonic drive［D］.Chengdu：Sichuan University，2005（in Chinese）.

［4］ KAYABASI O，ERZINCANLI F.Shape optimization of tooth profile of a flexspline for a harmonic drive by finite element modeling［J］.Materials&Design，2007，28（2）：441-447.

［5］秦磊，許立忠．機電集成靜電諧波傳動柔輪的應力分析［J］．機械傳動，2007，31（6）：13-15．QIN L，XU L Z.The stress analysis of the flexible of the electromechanical integrated electrostatic harmonic drive［J］.Mechanical Transmission，2007，31（6）：13-15（in Chinese）．

［6］王知行，徐克．杯型柔輪的應力分析與研究［J］．齒輪，1989，13（2）：23-26．WANG Z X，XU K.Stress analysis and study of cup-shaped flexspline［J］.Journal of Gear，1989，13（2）：23-26（inChinese）．

［7］韓建華，范祖堯．諧波傳動中杯形柔輪的應力分析 ［J］．齒輪，1987，11（5）：6-10. HAN J H，FAN Z R.Stress analysis and study of cup-shaped flexspline in harmonic drive［J］.Journal of Gear，1987，11（5）：6-10（in Chinese）.

［8］CHEN X X，LIU Y S，XING J Z，et al.A novel method based on mechanical analysis for the stretch of the neutral line of the flexspline cup of a harmonic drive［J］.Mechanism and Machine Theory，2014，76：1-19.

［9］CHEN Xiao-xia，LIN Shu-zhong，XING Jing-zhong.Modeling of flexspline and contact analyses of harmonic drive［J］. Key Engineering Materials，2010，419/420：597-600.

［10］付軍鋒，董海軍，沈允文.諧波齒輪傳動中柔輪應力的有限元分析［J］.中國機械工程，2007，18：2210-2214. FU J F，DONG H J，SHEN Y W.The finite element analysis of the stress of flexible in harmonic drive［J］.Chinese Mechanical Engineering，2007，18：2210-2214（in Chinese）.

［11］劉文芝，張乃人，張春林，等.諧波齒輪傳動中杯形柔輪的有限元計算與分析［J］.機械工程學報，2006，42（4）：52-57. LIU W Z，ZHANG N R，ZHANG C L，et al.The finite element calculation and analysis of the cup-shaped flexible in harmonic drive［J］.Journal of Mechanical Engineering，2006，42（4）：52-57（in Chinese）.

［12］吳偉國，張勇，梁風.具有嚙合齒面接觸對的諧波傳動有限元模型建立與分析［J］.機械傳動，2011，12：37-41. WU W G，ZHANG Y，LIANG F.The establishment and analysis of the finite element model of harmonic drive with the contact of meshing tooth surface［J］.Mechanical Transmission，2011，12：37-41（in Chinese）.

［13］DONG H，TING K L，WANG D.Kinematic fundamentals of planer harmonic drives［J］.ASME J Mesh Des，2011，133（1）：（011007-1）-（011007-8）.

［14］DONG H，WANG D，TING K L.Kinematic effect of the compliant cup in harmonic drives［J］.ASME J Mech Des，2011，133（1）：（051004-1）-（051004-7）.

［15］董惠敏，張曉青.基于實驗建模的諧波齒輪傳動柔輪的有限元分析研究［J］.機械傳動，200l，25（2）：16-19. DONG H M，ZHANG X Q.Finite element analysis of the flexible of harmonic with experimental model［J］.Mechanical Transmission，200l，25（2）：16-19（in Chinese）.

［16］陳曉霞.基于側隙評價的諧波齒輪圓弧齒廓設計及狀態(tài)仿真研究［D］.天津：天津工業(yè)大學，2012. CHEN X X.Design and simulation study of the arc tooth profile of harmonic based on backlash evaluation［D］.Tianjin：Tianjin Polytechnic University，2012（in Chinese）.

［17］機械設計手冊編委會.機械設計手冊單行本：輪系［M］.北京：機械工業(yè)出版社，2007：166-169. Mechanical Design Manual Editorial Board.Mechanical Design Handbook Offprint：Gear［M］.Beijing：Mechanical Industry Press，2007：166-169（in Chinese）.

Stress analysis of harmonic gear driven in assemble state with 3D solid finite element model

CHEN Xiao-xia1，2，CHENG Kai1
（1.School of Mechanical Engineering，Tianjin Polytechnic University，Tianjin 300387，China；2.Tianjin Key Laboratory of Modern Mechatronics Equipment Technology，Tianjin Polytechnic University，Tianjin 300387，China）

In order to describe the tooth profile and dedendum stress distribution of harmonic flexspline correctly，a 3D solid model method based on real tooth profile is put forward.Using involute tooth profile as an example，a parameterized solid model based on language APDL is set up，which includes important parameters such as involute tooth profile，tooth thickness and tooth root fillet.The stress of flexspline under the action of wave generator is obtained by defining contact between the wave generator and flexspline's inner surface.The stress is compared with the theoretical value based on equivalent ring theory.In order to analyze the stress distribution of flexspline in assemble state more fully，a 3D solid model of flexspline is built and solved by ANSYS.The stress is extracted by defining the path in the middle section corresponding to the section calculated by theoretical formula and compared with the results obtained by the planar gear ring and theoretical value.The results show that in the finite element model of planar gear ring，the mean value of the circumferential stress fluctuation is consistent with the theoretical value，which validates the validity of the finite element model method.In addition，the circumferential stress extracted from 3D solid model is larger than that from planar gear ring，and the maximum circumferential stress occurs at the top layer.

harmonic gear；flexspline；finite element analysis；3D solid model；assemble state

TH112.2

A

1671-024X（2016）05-0077-06

10.3969/j.issn.1671-024x.2016.05.014

2016-03-09

國家自然科學基金面上項目（51575390）；天津市應用基礎與前沿技術研究計劃一般項目（14JCYBJC19200）

陳曉霞（1969—），女，博士，副教授，主要研究方向為諧波齒輪傳動、數字化設計制造.E-mail：chenxx@tjpu.edu.cn