摘    要：基于ABAQUS-Python提出了一種橡膠-金屬襯套快速仿真技術。該技術將典型橡膠-金屬襯套結構參數化，并通過開發的獨立圖形用戶界面和ABAQUS腳本程序，實現自動前處理、仿真計算和后處理；讀取仿真結果文件中力、扭矩、位移和角度值，采用最小二乘法計算出多向靜剛度值，導出應力、應變等云圖；對比仿真與實測結果，誤差在10%以內，滿足工程化應用要求。此外，該方法進行一次仿真分析約需8～15 min,極大地提高了分析效率。

關鍵詞：參數化;橡膠-金屬襯套;仿真技術;

隨著汽車工業的迅速發展，越來越多的人們重視整車舒適性和操穩性。作為一種阻尼材料，橡膠具有良好的彈性特性和能量吸收能力，因此被廣泛應用于汽車減振領域，如底盤襯套、動力總成懸置和其他橡膠制品等。

底盤襯套類橡膠減振制品由芯軸、外套和橡膠組成，通過橡膠硫化過程，實現3者連接。常規橡膠襯套可通過調節各組件結構、尺寸和橡膠硬度，實現在x、y、z軸平移和偏轉的性能要求。由于車型和車輛系統各部位的差異很大，對橡膠減振制品的性能要求也不同，因此需根據不同的性能要求選擇合適的結構、尺寸和橡膠材料，并進行優化設計。

在設計橡膠襯套類制品時，常采用試驗法和數值模擬分析法來確定相關參數。數值模擬分析法因具有直觀、快速、成本低等優點而被廣泛應用，但數值模擬分析法中的建模、前處理和后處理等過程繁瑣，費時費力。因此，研究參數化建模、自動前處理和后處理的方法，對橡膠襯套類制品的快速開發和優化有顯著的實際工程意義。

為了提高橡膠-金屬件數值模擬分析的效率，劉志國等[1]通過模型重建與專家分析經驗封裝構建軌道車輛金屬橡膠件有限元模型參數化建模方法，實現了金屬橡膠件分析模型的參數化與仿真分析自動化。周煒等[2]利用剛度疊加原理和參數有限元分析，建立了橡膠襯套幾何參數與剛度的關系式；針對剛度匹配的設計要求，通過設定不同的目標函數，得到了橡膠襯套結構的最優解。宋穎等[3]基于ABAQUS二次開發功能，開發了圓柱形橡膠堆的插件程序，實現了圓形橡膠堆參數化建模和后處理的自動化。

綜上所述，國內眾多學者已對橡膠金屬襯套件的參數化建模和快速仿真進行了研究，并且利用現有的有限元軟件進行了二次開發，提高了開發效率。但這些研究都無法實現輸入參數后快速得到結果的目的。鑒于此，本文基于ABAQUS-Python提出了一種橡膠-金屬襯套件快速仿真分析技術。該方法構建了獨立的交互界面，無需啟動ABAQUS有限元軟件，用戶只需輸入產品結構參數，即可快速得到靜剛度值并輸出應力、應變等結果云圖。本文的研究成果為提高橡膠-金屬襯套件開發效率打下了堅實的基礎，具有操作簡單、快速輸出結果等優點。

1 橡膠-金屬襯套結構參數化建模

1.1 尺寸結構參數提取

將橡膠-金屬襯套結構進行參數化，包括襯套的尺寸、形狀等參數，如圖1所示。具體參數詳見表1。

圖1 橡膠-金屬襯套結構和尺寸

表1 橡膠-金屬件尺寸和材料參數

芯軸的外側、外套的內外側通常為圓柱形或旋轉面。芯軸和外套通常采用金屬或塑料制成，剛度明顯高于橡膠。為了簡化分析，提高效率，在前處理過程中省略芯軸和外套，并將橡膠內表面耦合到一個點。此外，為了實現襯套結構中外套的縮頸過程，將橡膠外表面定義為柱面坐標系。

1.2 基于ABAQUS-Python構建程序

1.2.1 參數化建模

本文研究的橡膠-金屬襯套結構為回轉體結構，上下對稱，因此只需繪制1/4結構，并通過鏡像的方法得到完整模型。圖2展示了繪制的1/4結構。

圖2 橡膠-金屬襯套1/4結構

結合表1中的數據，運用幾何計算的方法，易得圖2中6個點和其他各點坐標值。本文編寫的程序可得每條特征線段和圓弧，其中部分直線和圓弧代碼為：

s.Line(point1=(X2, Y2), point2=(X4, Y4)) #畫出直線

s.ArcByCenterEnds(center=(X3, Y3), point1=(X2, Y2), point2=(X1, Y1), direction=CLOCKWISE) #畫出圓弧

為保證生成的橡膠型面沒有線條交叉和尖點，圖2中6個點坐標必須滿足式(1)所示約束條件。

1.2.2 定義橡膠本構

超彈理論認為：橡膠材料具有與應變率無關的非線性彈性特性，隨著外力的消失，材料可以恢復到原始狀態，橡膠材料應力應變關系可以用單位體積應變能來度量。Yeoh[4]通過分析應變能函數，提出了三次應變能函數，如式(2)所示，稱為Yeoh本構方程。

W=C10 (I1-3)+C20 (I2-3)2+C30 (I1-3)3 (2)

式中：W為應變能函數，C10、C20、C30為材料常數；I1為Cauchy-Green變形張量的第一基本不變量。

Yeoh本構方程得到的應變能與試驗結果較吻合，可由簡單的單軸拉伸試驗擬合出材料參數，進而用來預測其他的力學行為，常用于橡膠制品結構的大變形分析[5,6]。本文采用Yeoh本構，參數為實測結果經擬合得出。部分代碼如下：

mymaterial = myModel.Material(name=′Rubber_′+str(hardness) +′ShA′)

mymaterial.Density(dependencies=1, table=((1.2e-09, 0.0),))

mymaterial.Hyperelastic(materialType=ISOTROPIC, testData=OFF, type=YEOH, volumetricResponse=VOLUMETRIC_DATA, table=((c10, c20, c30, d1, d2, d3),))

mymaterial.Expansion(table=((expantemp, 0.0),), dependencies=1)

1.2.3 定義邊界和輸出類型

為了保證分析精度，同時提高分析效率，本文省略了橡膠襯套的內、外套等不必要細節，橡膠內圓柱面采用rbe2方式耦合為一點，并將橡膠外圓柱面定義為柱坐標系。

分析過程主要分為兩步：第1步，耦合點固定，外圓柱面第一坐標值的變化量為(L-F)/2;第2步，外圓柱面固定，耦合點于徑向、軸向、偏擺和扭轉方向分別加載0.5 mm、0.5 mm、3°和3°。

輸出類型定義為常見的應力、應變、位移、扭矩、支反力等，部分代碼如下：

myModel.FieldOutputRequest(name=′F-Output-1′, createStepName=′cool_ass′, variables=(′S′,′LE′,′U′,′UR′,′RF′,′RM′,′CF′,′CSTRESS′,′TEMP′))

myModel.HistoryOutputRequest(name=′H-Output-1′, createStepName=′cool_ass′, variables=(′U1′,′U2′,′U3′,′UR1′,′UR2′,′UR3′,′RF1′,′RF2′,′RF3′,′RM1′,′RM2′,′RM3′,′CF1′,′CF2′,′CF3′,′CM1′,′CM2′,′CM3′), region=regionDef, sectionPoints=DEFAULT, rebar=EXCLUDE)

1.2.4 網格劃分和提交計算

為提高計算的速度和精度，本文采用六面體網格，并且選擇了C3D8H類型的網格。利用submit語句完成計算的提交，waitForCompletion()語句可實現多工況的順序計算。橡膠-金屬制品網格模型如圖3所示。網格劃分和提交計算的部分代碼如下：

partInstances = (myAssembly.instances[′rubberInstance′],)

myAssembly.seedPartInstance(regions=partInstances, size=1.0, deviationFactor=0.1, minSizeFactor=0.1)

myAssembly.generateMesh(regions=partInstances)

elemType = mesh.ElemType(elemCode=C3D8H, elemLibrary=STANDARD)

myJob = mdb.Job(name=jobName, model=modelName, type=ANALYSIS, numCpus=numCores)

myJob.submit()

myJob.waitForCompletion()

圖3 橡膠-金屬襯套網格模型

1.2.5 后處理

計算完成后，可讀取相應的ODB文件提取位移或旋轉角度、支反力或扭矩等信息，輸出應力、應變云圖；利用最小二乘法計算各向靜剛度值，并繪制曲線。以徑向計算為例，圖4示出了橡膠最大主應力結果，圖5示出了橡膠最大主應變結果，圖6示出了位移與支反力曲線。靜剛度輸出涉及部分代碼如下：

def LSM(data1,data2): #最小二乘法

lenData = len(data1)

averageX = sum(data1)/lenData; averageY= sum(data2)/lenData

upperData = 0; lowerData = 0

for i in range(lenData):

upperData += (data1[i]-averageX)*(data2[i]-averageY)

lowerData += (data1[i]-averageX)**2

stiffness = upperData/lowerData

inaccuracy = averageY-averageX*b

return stiffness, inaccuracy

tempDataX=myOdb.steps[loadCaseName].historyRegions[′Node ASSEMBLY.1′].historyOutputs[′U1′]

for i, j in tempDataX.data: dataX.append(j)

tempDataY=myOdb.steps[loadCaseName].historyRegions[′Node ASSEMBLY.1′].historyOutputs[′RF1′]

for i, j in tempDataY.data: dataY.append(j)

stiffnessAZ, inaccuracyDataAZ=LSM(dataX, dataY)

圖4 最大主應力

圖5 最大主應變

圖6 位移和支反力曲線

2 開發獨立圖形用戶界面

圖形用戶界面是一種將人、機、環境3者作為一個系統進行總體設計，強調人性化、操作方便等特點的人機系統工程。它對使用者的認知要求較低，更適合普通使用者操作使用。

圖7 獨立圖形用戶界面

本文第1節中開發的程序可以通過ABAQUS中‘Run Script’或 ‘Plug-ins’運行，但需要先運行ABAQUS。為了進一步節約時間、降低使用難度，本文利用python中tkinter庫，利用選值框、畫布和邊框等指令開發獨立于ABAQUS的圖形用戶界面，如圖7所示，其中右側的橡膠截面可隨輸入參數的變化同步更新。

當點擊‘Submit’時，程序啟動運行run.bat文件，調用bat文件的代碼如下：

thisPath =os.path.dirname(os.path.abspath(__file__))

file_path = thisPath +′\\′+′run.bat′

Popen(file_path, creationflags=subprocess.CREATE_NEW_CONSOLE)

文件run.bat用于在后臺運行第1節編寫的程序(normalFEA.py),run.bat文件代碼如下：

abaquscae noGUI=normalFEA.py

3 試驗對比和分析驗證

橡膠襯套剛度測試采用INOVA雙軸測試機，型號為FEH 25-7-400 v1,產品測試見圖8。本文選取3種不同尺寸的橡膠-金屬襯套制品，產品的結構尺寸和橡膠參數詳見表2。

圖8 產品測試

仿真和實測結果以及誤差值見表3,其中誤差值取(仿真值-實測值)與仿真值的比值。通過表3可知，仿真值與實測值間存在一定的誤差。這是部分由于實際生產和仿真分析中存在誤差導致的，比如實際生產過程中膠料硬度會存在波動，而仿真分析中采用的橡膠硬度是名義值；仿真分析中對結構的簡化和網格參數的選取會影響計算的準確度。誤差值控制在10%以內，滿足工程化應用的要求。經過大量案例仿真的用時統計，本文涉及的橡膠-金屬襯套類結構，絕大多數產品的仿真分析可在8～15 min內完成，極大地提高了分析效率。

表2 產品參數

表3 仿真和實測值

4 結論

本文基于ABAQUS-Python接口，進行ABAQUS二次開發；開發獨立于ABAQUS的圖形用戶界面，得出如下主要結論。

(1) 該方法和過程可實現橡膠-金屬襯套制品的有限元快速分析。前處理、計算、后處理等工作可自動完成，仿真分析可在8～15 min內完成，極大地提高了分析效率。

(2) 通過與實測結果分析對比，結果誤差小于10%,滿足工程化應用的要求。

(3) 操作簡單、高效，使用者無需具有專業的有限元分析知識即可進行仿真分析，極大地降低了使用門檻。

參考文獻

[1] 劉志國，王炳杰.軌道車輛金屬橡膠件參數化建模及仿真分析方法研究[J].鐵道車輛，2017,55(12):8-15.

[2] 周煒，黃友劍，李建林.基于參數化有限元的橡膠襯套結構優化設計[J].特種橡膠制品，2012,33(4):50-54.

[3] 宋穎，榮繼剛，馬秋成，等.基于ABAQUS的圓柱形橡膠堆參數化建模的二次開發[J].信息技術，2022,52(1):146-148.

[4] YEOH OH.Some Forms of the Strain Energy Function for Rubber[J].Rubber Chemistry and Technology,1993,66(5):754-771.

[5] 李曉芳，楊曉翔.橡膠材料的超彈性本構模型[J].彈性體，2005,15 (1):50-58.

[6] 朱艷峰，劉鋒，黃小清，等.橡膠材料的本構模型[J].橡膠工業，2006,53(2):119-125

文章來源特種橡膠制品. 2023,44(05)