ABAQUS二次開發(fā)在直齒輪彎曲應(yīng)力分析中的應(yīng)用

王玉玲，張榮闖，李明

（1.沈陽(yáng)城市建設(shè)學(xué)院機(jī)械工程學(xué)院，遼寧 沈陽(yáng) 110167；2.東北大學(xué)秦皇島分校控制工程學(xué)院，河北 秦皇島 066004）

齒輪作為機(jī)械產(chǎn)品傳動(dòng)部件應(yīng)用中最為廣泛的零件，對(duì)其進(jìn)行彎曲應(yīng)力分析是保證使用壽命和傳動(dòng)能力的前提。目前，齒輪彎曲應(yīng)力計(jì)算的各類計(jì)算標(biāo)準(zhǔn)中，大多數(shù)基于Lewis公式，該公式基于材料力學(xué)等強(qiáng)度懸臂梁假設(shè)，沒有考慮齒根截面的突變，同時(shí)忽略了輪齒徑向載荷對(duì)應(yīng)力的影響，導(dǎo)致不能準(zhǔn)確地反映齒輪的彎曲應(yīng)力［1-3］。有限元分析法的應(yīng)用，能夠準(zhǔn)確地描述輪齒的幾何形狀、邊界條件以及載荷狀況，使齒輪彎曲應(yīng)力計(jì)算更加便捷準(zhǔn)確。然而采用CAE軟件進(jìn)行齒輪彎曲應(yīng)力分析，圖形界面操作過(guò)程繁冗復(fù)雜，仿真分析效率較低［4-7］。有限元分析軟件ABAQUS采用基于特征的參數(shù)化建模方式，具有基于Python腳本編程的二次開發(fā)功能，其高效性和便捷性已得到驗(yàn)證［8-9］。本文采用Python語(yǔ)言編寫內(nèi)核腳本程序，從而控制ABAQUS內(nèi)核操作，實(shí)現(xiàn)齒輪彎曲應(yīng)力有限元分析的參數(shù)化和自動(dòng)化。

2 輪齒幾何建模

采用三齒模型進(jìn)行齒輪彎曲應(yīng)力分析。齒廓幾何模型采用文獻(xiàn)［10］的圓弧過(guò)渡曲線齒廓模型，如圖1所示，坐標(biāo)系OG-XGYG為齒廓坐標(biāo)系，原點(diǎn)OG位于齒輪中心，坐標(biāo)軸OGYG與齒廓中心線重合，rf、rb、r以及ra分別為齒根圓、基圓、分度圓和齒頂圓的半徑，過(guò)渡曲線采用圓弧曲線，其分別與漸開線、齒根圓弧相切于B點(diǎn)和D點(diǎn)，圓心位于E點(diǎn)。

圖1 齒廓幾何模型Fig.1 Tooth profile geometry model

齒廓計(jì)算的各關(guān)鍵點(diǎn)坐標(biāo)計(jì)算公式如下：

式中：ΩS為1/2基圓齒厚對(duì)應(yīng)的圓心角，ΩS=π/2z+θn，θn=tanαn-αn，αn為分度圓壓力角，z為齒數(shù)；SR為輪緣厚度，SR=htmB，mB為輪緣齒高比，SR的大小直接影響齒根彎曲應(yīng)力，根據(jù)標(biāo)準(zhǔn)GB/T/ISO 6336-3∶2019，mB取值為1.2，此時(shí)可忽略輪緣厚度對(duì)齒根彎曲應(yīng)力的影響。

式中：αk、θk分別為K點(diǎn)對(duì)應(yīng)的壓力角和展角，θk=tanαk-αk，αk的取值為0≤αk≤arccos

3 有限元模型的建立

在低速重載齒輪傳動(dòng)中，需要按照彎曲強(qiáng)度計(jì)算齒輪的承載能力，即分析齒輪危險(xiǎn)截面處的最大彎曲應(yīng)力。假設(shè)齒輪傳動(dòng)為勻速傳動(dòng)，忽略輪齒誤差。z1齒輪的精度等級(jí)為6級(jí)，計(jì)算齒根基本應(yīng)力值σF0為44 MPa，應(yīng)力修正系數(shù)YS為1.97；轉(zhuǎn)速為200 r/min，轉(zhuǎn)矩為120 kN·m。齒輪的相關(guān)幾何參數(shù)見表1，z1齒輪的材料特性見表2。

表1 齒輪幾何參數(shù)Tab.1 Gear geometric parameters

表2 z1齒輪材料特性Tab.2 z1 gear material

3.1 邊界條件與載荷

常用直齒輪傳動(dòng)重合度系數(shù)一般在1～2之間，傳動(dòng)過(guò)程中會(huì)出現(xiàn)單齒、雙齒交替嚙合現(xiàn)象；當(dāng)輪齒在單對(duì)齒嚙合區(qū)外界點(diǎn)嚙合時(shí)，只有一對(duì)輪齒承載，此時(shí)齒根所受的彎矩最大。因此，齒輪的齒根彎曲應(yīng)力應(yīng)該按載荷作用于單對(duì)齒嚙合區(qū)外界點(diǎn)來(lái)計(jì)算。如圖2所示，B1為單對(duì)齒嚙合區(qū)外界點(diǎn)，通過(guò)式（10）可以求取坐標(biāo)系OG-XGYG下的載荷作用點(diǎn)（x0，y0）。

圖2 單對(duì)齒嚙合區(qū)外界點(diǎn)計(jì)算Fig.2 Determination of highest point sign tooth contact

αFen載荷作用角用來(lái)確定集中力載荷Fn的方向：

假設(shè)齒輪運(yùn)動(dòng)傳遞扭矩為T，單對(duì)齒嚙合區(qū)外界點(diǎn)集中力載荷Fn為

位移約束施加在左右輪緣和底部輪緣，約束全部位移。載荷施加與邊界條件定義如圖3所示。

圖3 載荷施加與邊界條件定義Fig.3 Load application and boundary condition definition

3.2 網(wǎng)格劃分

網(wǎng)格劃分是否合理，直接影響著有限元分析計(jì)算的精度和效率；為了保證齒輪網(wǎng)格劃分的質(zhì)量以及參數(shù)化實(shí)現(xiàn)，采用分塊劃分的策略；根據(jù)齒廓的幾何形狀、載荷作用點(diǎn)的位置，將單個(gè)輪齒切割為8個(gè)區(qū)域塊，如圖4所示。三齒模型的分塊結(jié)果如圖5所示。

圖4 單齒模型區(qū)域劃分Fig.4 Region division for single tooth

圖5 三齒模型區(qū)域劃分Fig.5 Region division for three teeth

單元形狀四邊形，單元類型為平面應(yīng)變單元CPE4R，采用結(jié)構(gòu)網(wǎng)格劃分技術(shù)；網(wǎng)格劃分時(shí)考慮齒的工作部分應(yīng)力集中和計(jì)算效率，嚙合齒廓和過(guò)渡曲線區(qū)域采用比較密集的網(wǎng)格，輪齒輪轂基體采用稀疏網(wǎng)格；網(wǎng)格疏密程度是通過(guò)控制區(qū)域塊邊界控制線的網(wǎng)格劃分種子數(shù)量來(lái)實(shí)現(xiàn)的。

網(wǎng)格劃分后的有限元模型如圖6所示，整個(gè)模型共劃分單元25 200個(gè)，節(jié)點(diǎn)25 731個(gè)。

圖6 網(wǎng)格劃分Fig.6 Mesh devision

4 仿真模型的代碼實(shí)現(xiàn)

ABAQUS采用Python作為其腳本命令開發(fā)語(yǔ)言，在Python原有編程對(duì)象的基礎(chǔ)上，融合了ABAQUS特有的Session、Mdb和Odb 3個(gè)對(duì)象。Session對(duì)象負(fù)責(zé)視圖的控制與定義；Mdb對(duì)象負(fù)責(zé)齒輪的幾何建模、材料定義、網(wǎng)格劃分以及邊界條件定義等功能；Odb對(duì)象負(fù)責(zé)仿真結(jié)果的輸出。

4.1 幾何建模

由圖1可知，齒頂圓、過(guò)渡曲線、齒根圓和底部輪緣均為圓弧曲線，可以通過(guò)圓心點(diǎn)和兩端點(diǎn)表示，例如表示為圓心E、起點(diǎn)B和終點(diǎn)D的圓弧曲線。輪緣FG為直線，可以通過(guò)起點(diǎn)F和終點(diǎn)G來(lái)表示。漸開線可根據(jù)式（9），表示為多個(gè)離散幾何點(diǎn)擬合而成的樣條曲線。對(duì)應(yīng)于三齒齒廓可以通過(guò)對(duì)應(yīng)幾何曲線的坐標(biāo)變換來(lái)求取。上述齒廓的幾何模型形狀取決于齒輪模數(shù)、齒數(shù)以及分度圓壓力角，因此可實(shí)現(xiàn)參數(shù)化建模。幾何建模的關(guān)鍵腳本代碼如下：

s=mdb.models['Model-1'].Constrained Sketch(name='__profile__',sheetSize=200.0)

s.Line(point1=(),point2=())#兩點(diǎn)直線構(gòu)造側(cè)輪緣;

s.ArcByCenterEnds(center=(),point1=(),point2=(),direction=CLOCKWISE)#中心點(diǎn)和兩端點(diǎn)構(gòu)造過(guò)渡曲線、齒根圓、齒頂圓和底部輪緣;

s.Spline()#樣條曲線構(gòu)造漸開線;

p=mdb.models['Model-1'].parts['Part-1']

f=p.faces

picked Faces=f.find At(coordinates=(0.0,r,0.0))

p.PartitionFaceBySketch(faces=pickedFaces,sketch=s)#輪齒區(qū)域劃分主要在草圖中通過(guò)直線分割面功能來(lái)實(shí)現(xiàn)。

4.2 材料定義

mdb.models['Model-1'].Material(name='Material-1')

mdb.models['Model-1'].materials['Material-1'].Elastic(table=((210000.0,0.3),))#定義彈性模量和泊松比;

mdb.models['Model-1'].Homogeneous Solid Section(name='Section-1',material='Material-1',thickness=None)#創(chuàng)建截面;

region=p.Set(faces=f,name='Set-1')

p.SectionAssignment(region=region,sectionName='Section-1',offset=0.0,offsetType=MIDDLE_SURFACE,offsetField='',thicknessAssignment=FROM_SECTION)#分配截面。

4.3 施加約束和載荷

mdb.models['Model-1'].DisplacementBC(name='BC-1',createStepName='Initial',region=region,u1=SET,u2=SET,ur3=UNSET,amplitude=UNSET,

distributionType=UNIFORM,fieldName='',localCsys=None)#邊界條件定義,施加約束;

KLM=a.DatumCsysByThreePoints(name='Datum csys-20',coordSysType=CARTESIAN,origin=(),point1=(),point2=())#載荷為集中力,集中力作用于單對(duì)齒嚙合區(qū)外界點(diǎn),方向垂直于漸開線,需要定義局部坐標(biāo)系確定集中力方向,三點(diǎn)定義局部坐標(biāo)系;

datum=mdb.models['Model-1'].rootAssembly.datums[KLM.id]

v1=a.instances['Part-1-1'].vertices verts1=v1.find-At(((),))#載荷作用點(diǎn);

region9=a.Set(vertices=verts1,name='Set-9')

mdb.models['Model-1'].ConcentratedForce(name='Load-1',createStepName='applyload',region=region9,cf1=-1800.0,distributionType=UNIFORM,field='',localCsys=datum)#集中力施加。

4.4 網(wǎng)格劃分

datumP=p.Datum Point By Coordinate(coords=())

d1=p.datums

p.Partition Edge By Point(edge=e1.find At(coordinates=)),point=d1[datumP.id])#網(wǎng)格的疏密程度通過(guò)指定邊的種子數(shù)來(lái)確定。方便于不同區(qū)域網(wǎng)格的劃分,需要將邊以通過(guò)點(diǎn)的方式進(jìn)行分割處理;

pickedEdges=e.findAt(((),))

p.seedEdgeByNumber(edges=pickedEdges,number=10,constraint=FIXED)#設(shè)定邊的種子數(shù)量;

pickedRegions=p.faces

p.setMeshControls(regions=pickedRegions,elemShape=QUAD,technique=STRUCTURED)

elemType1=mesh.ElemType(elemCode=CPE4R,elemLibrary=STANDARD,

secondOrderAccuracy=OFF,hourglassControl=DE-FAULT,

distortionControl=DEFAULT)

elemType2=mesh.ElemType(elemCode=CPE3,elemLibrary=STANDARD)#設(shè)定單元類型;

p=mdb.models['Model-1'].parts['Part-1']

f1=p.faces

pickedRegions1=(f1,)

p.setElementType(regions=pickedRegions1,elemTypes=(elemType1,elemType2))

p.generateMesh()#網(wǎng)格劃分。

5 仿真結(jié)果分析與對(duì)比

仿真結(jié)果如圖7所示，左側(cè)過(guò)渡曲線彎曲應(yīng)力最大值為47.000 0 MPa，右側(cè)過(guò)渡曲線彎曲應(yīng)力最大值為41.559 7 MPa，鑒于載荷作用點(diǎn)為單對(duì)齒嚙合區(qū)外界點(diǎn)應(yīng)力修正系數(shù)Ys的影響，齒輪計(jì)算彎曲應(yīng)力為88.65 MPa，有限元仿真與計(jì)算值相差較大，其主要原因是國(guó)際標(biāo)準(zhǔn)計(jì)算結(jié)果明顯過(guò)于保守，同時(shí)采用圓弧作為過(guò)渡曲線，其抗彎和抗拉能力要優(yōu)于外擺線的過(guò)渡曲線。

6 結(jié)語(yǔ)

本文準(zhǔn)確描述圓弧過(guò)渡曲線的直齒輪齒廓數(shù)學(xué)模型，重點(diǎn)描述齒輪彎曲應(yīng)力有限元建模過(guò)程中載荷施加作用點(diǎn)位置的計(jì)算方法以及邊界條件的施加方式。根據(jù)網(wǎng)格劃分的疏密要求，研究了齒廓幾何模型的區(qū)域劃分方法和網(wǎng)格種子分布方式。給出了齒輪彎曲應(yīng)力有限元仿真模型基于Python的ABAQUS二次開發(fā)編程代碼，實(shí)現(xiàn)了應(yīng)力分析過(guò)程的參數(shù)化和自動(dòng)化。應(yīng)用基于Python的ABAQUS二次開發(fā)直齒輪彎曲應(yīng)力分析，為齒輪傳動(dòng)過(guò)程其他物理特性的高效便捷有限元計(jì)算（例如接觸分析）提供了借鑒。