吳 躍，王廷劍，，3，孫厚濤，齊智慧，張永釗，王燕霜，宣兆成

（1.天津職業(yè)技術(shù)師范大學機械工程學院，天津 300222；2.哈爾濱工業(yè)大學機電工程學院，哈爾濱 150001；3.中國航發(fā)哈爾濱軸承有限公司，哈爾濱 150001）

隨著材料技術(shù)和表面工程技術(shù)的發(fā)展，固體潤滑涂層技術(shù)被廣泛應(yīng)用于齒輪和軸承等關(guān)鍵基礎(chǔ)件的工作表面，以提高其減摩、耐磨、抗疲勞和抗膠合等綜合摩擦學性能，特別是工作在高真空、強輻射和高低溫交替的空間苛刻環(huán)境的摩擦部件[1-3]。除了涂層材料自身的機械力學性能及其與基體的界面結(jié)合強度之外，涂層在接觸載荷作用下的力學行為也是決定涂層實際應(yīng)用效果的核心因素，因此涂層體系的接觸力學行為是涂層優(yōu)化設(shè)計不可忽略的問題。

許多研究學者采用ANSYS、ABAQUAS 等商用有限元軟件分析齒輪、軸承等接觸界面的力學行為，解釋試驗結(jié)果，揭示涂層的失效機制[4-5]。應(yīng)用有限元軟件分析接觸問題，其計算結(jié)果的精度和準確性依賴于模型網(wǎng)格單元劃分質(zhì)量以及接觸剛度的設(shè)置等，這對于普通工程技術(shù)人員而言非常不便。積分變換法是彈性力學應(yīng)力應(yīng)變分析常用的方法，文獻[6-9]應(yīng)用積分變換的方法獲得了單層及多層涂層體系半平面在表面線分布機械載荷作用下彈性場的頻域解，建立了涂層體系的線接觸模型；文獻[10-14]應(yīng)用積分變換的方法獲得了單層和多層涂層體系半空間在表面機械載荷作用下彈性場的頻域解，建立了涂層體系的三維點接觸模型，并基于與快速傅里葉變換相關(guān)的多種算法開發(fā)了相應(yīng)的求解分析程序，該程序具有計算速度快、計算精度高等優(yōu)點，但其缺乏軟件用戶界面，以致接觸參數(shù)輸入和計算結(jié)果查看不方便。在工程領(lǐng)域，通過Vb、Fortran 和Matlab 等計算機語言或應(yīng)用軟件間的混合編程，為各種工程問題的專用計算分析程序開發(fā)軟件界面，方便工程技術(shù)人員使用是一種發(fā)展趨勢。劉鵬等[15]應(yīng)用Vb 中的Shell 功能函數(shù)直接調(diào)用Fortran 外殼程序，開發(fā)了結(jié)構(gòu)動力分析有限元計算程序系統(tǒng)。郭翔偉等[16]通過動態(tài)鏈接庫.dll 文件實現(xiàn)Vb對Fortran 混合編程，開發(fā)了混凝土壩溫度應(yīng)力仿真分析軟件。萬庭輝等[17]利用ActiveX 技術(shù)實現(xiàn)Vb 與Matlab 的混合編程，開發(fā)了TOUGH+HYDRATE 數(shù)值模擬軟件。馬培勇等[18]利用COM 組件生成.dll 文件來實現(xiàn)Vb 對Matlab 混編，開發(fā)了圓柱形彈簧變形量的蒙特卡洛模擬分析軟件。大量相關(guān)文獻表明，Vb、Fortran 和Matlab 等計算機語言或應(yīng)用軟件間的混合編程可以綜合利用Fortran 計算速度快、Vb 界面開發(fā)功能強大和Matlab 軟件繪圖功能強大的優(yōu)勢。

基于前期研究開發(fā)的單層涂層體系接觸分析的Fortran 程序代碼[13-14]，本文通過Vb、Fortran 和Matlab的混合編程，開發(fā)了單層涂層體系三維點接觸分析軟件——CONTACT3D VERSION 1.0.exe，實現(xiàn)了接觸參數(shù)輸入界面化和接觸分析結(jié)果界面可視化，并對軟件進行演示與計算結(jié)果對比驗證。

1 基于Vb 的軟件界面設(shè)計

運行Vb 6.0 軟件，創(chuàng)建一個新工程，將其窗口Form1的Caption 屬性設(shè)置為CONTACT Version 1.0。然后在該窗口添加3 個Frame 控件，并將3 個Frame 控件的Caption 分別設(shè)置為數(shù)值計算、接觸壓力分布和馮米斯等效應(yīng)力分布，Vb 開發(fā)的軟件界面如圖1 所示。

圖1 Vb 開發(fā)的軟件界面

在Caption 屬性為數(shù)值計算的Frame 控件下部創(chuàng)建一列TextBox 文本框，作為單層涂層體系與彈性球體間三維點接觸的接觸載荷、球體與涂層體系彈性參數(shù)、涂層厚度、幾何尺寸、摩擦系數(shù)及x、y、z 3 個方向劃分單元個數(shù)等的輸入，對應(yīng)的在其左側(cè)有一列Lable 控件，用于注釋對應(yīng)的文本框輸入?yún)?shù)的物理量和單位等內(nèi)容。在該Frame 控件上部添加一個Image 控件，并通過其Picture 屬性添加單層涂層體系半空間與彈性球體之間的三維點接觸示意圖，使界面輸入?yún)?shù)與其所代表的物理參量形成對照。在該Frame 控件最下方插入一個名稱為Command1 的命令按鈕，其Caption 屬性為開始計算，點擊該命令按鈕，其后臺程序?qū)⒄{(diào)用Fortran 編譯的單層涂層體系三維點接觸的Forlib.dll 動態(tài)數(shù)據(jù)庫里的DRYCONTACT過程稿完成接觸計算，并把接觸壓力和次表層馮米斯等效應(yīng)力的計算結(jié)果傳遞給Vb 的矩陣變量，同時以.txt 文件對包括接觸壓力和次表層馮米斯等效應(yīng)力在內(nèi)的數(shù)據(jù)進行保存，計算結(jié)果數(shù)據(jù).txt 文件存儲如圖2所示。圖中outputpressure.txt 為接觸壓力數(shù)據(jù)保存文件，VM.txt 為次表層等效馮米斯等效應(yīng)力的保存文件，其他.txt 文件為次表層接觸應(yīng)力6 個分量、兩接觸表面法向間距、接觸表面法向變形以及輸入的接觸工況參數(shù)等數(shù)據(jù)的保存文件。

圖2 計算結(jié)果數(shù)據(jù).txt 文件存儲

在Caption 屬性分別為接觸壓力分布和馮米斯等效應(yīng)力分布的Frame 控件內(nèi)分別添加一個Image 控件，其Name 屬性軟件默認名為Image1 和Image2，并設(shè)置好其大小。在2 個Image 控件下分別添加命令按鈕控件Command2 和Command3，命令按鈕的Caption屬性分別設(shè)置為接觸壓力和馮米斯等效應(yīng)力，通過點擊命令按鈕，其后臺程序?qū)⑼ㄟ^調(diào)用Matlab 的繪圖函數(shù)繪制接觸壓力分布和次表層馮米斯等效應(yīng)力梯度分布，實現(xiàn)計算結(jié)果的界面可視化。

2 Vb 和Fortran 的混合編程

2.1 動態(tài)鏈接庫.dll文件的編譯

Vb 與Fortran 的混合編程通常是通過Vb 調(diào)用Fortran 編譯的可執(zhí)行文件.exe 或者動態(tài)鏈接庫.dll 文件這2 種形式來實現(xiàn)[19]。動態(tài)鏈接庫文件是可供Vb、VC 等其他程序或者Fortran 自身調(diào)用的一些過程或函數(shù)組成的可執(zhí)行命令模塊，可以包含一個函數(shù)或者一個過程，也可以是同時包含多個函數(shù)或者多個過程，但其本身不能單獨運行。本文采用Vb 程序調(diào)用Fortran 編譯的動態(tài)鏈接庫文件的方法來實現(xiàn)單層涂層三維點接觸的分析計算，在前期研究開發(fā)的單層涂層體系接觸分析的Fortran 程序代碼[13-14]的基礎(chǔ)上，經(jīng)過簡要的形式改編后編譯成可供Vb 程序調(diào)用的動態(tài)鏈接庫文件。在Microsoft Visual Studio 2010 開發(fā)環(huán)境下，編譯動態(tài)鏈接庫文件的基本步驟如下。

步驟1打開Visual Studio Studio 2010，點擊菜單欄上的“文件”選項，選擇“新建”，然后選擇“項目”，彈出“新建項目”對話框，在對話框左側(cè)的模板列表中點擊Intel（R）VisualFortran 下的“Library”選項，然后在左側(cè)的Fortran 程序類型選項欄中選擇“Dynamic-link Library”，最后在下方輸入項目名和保存路徑，本文項目名為“Forlib”，完成后點擊確定按鈕。

步驟2在新出現(xiàn)的界面的右側(cè)，選擇sources files 選項，點擊“添加”按鈕，再點擊“新建項”，彈出一個新界面，在該界面里選擇“Fortran free-form file”選項，新建一個.f90 文件，本文將該Fortran 程序文件命名為“DRYCONTACT.f90”，然后點擊添加按鈕。

步驟3在開發(fā)環(huán)境窗口打開DRYCONTACT.f90文件，在這一步首先將涂層體系三維點接觸分析Fortran 程序代碼輸入到該文件，然后根據(jù)編譯.dll 動態(tài)數(shù)據(jù)庫的編程規(guī)則改編成函數(shù)或過程的形式，在改編過程中還需特別注意Vb 調(diào)用Fortran 編譯的.dll 文件的函數(shù)或過程的有關(guān)聲明，包括函數(shù)、過程或者變量的命名聲明、參數(shù)數(shù)據(jù)類型及參數(shù)傳遞方式等方面的約定。本文改編為過程的形式，改編完成后，在界面右側(cè)的sources files 選項下添加DRYCONTACT.f90 文件調(diào)用的各個子過程的f90 文件，包括G_INITIAL.f90、NCS.f90、ICONVCIJ.f90 和STRESSSOLVER.f90等，然后在菜單欄中選擇“生成”選項，選擇“生成解決方案”，最后選擇保存。

步驟4在該項目存儲路徑下的.debug 文件夾可以找到相應(yīng)的動態(tài)鏈接庫文件，本文動態(tài)鏈接庫文件為Forlib.dll。

在步驟3 中，DRYCONTACT.f90 文件的部分程序代碼見附錄1。程序中“! MS$ATTRIBUTES”是Microsoft Visual Studio 2010 Fortran 的元命令，也可寫為“! EDS$ATTRIBUTES”或“$ATTRIBUTES”，若寫成“$ATTRIBUTES”，則在Fortran 程序代碼中$必須位于首列。元命令主要用于聲明微軟擴展屬性，DLLEXPORT、ALIAS 和REFERENCE 都是微軟擴展屬性，DLLEXPORT 用于聲明過程或者函數(shù)可以被其他程序調(diào)用；ALIAS 用于為過程或者函數(shù)指定調(diào)用方的別名，在這由于調(diào)用方的過程名與Fortran 子過程的名字一致，因此也可以省略不寫；REFRENCE 用于聲明參數(shù)傳遞為按址（引用）傳遞，按址傳遞是Vb 與Fortran 混編默認的參數(shù)傳遞方式，通常也可以省略不寫。另一種常用的參數(shù)傳遞方式為按值傳遞，用VALUE 來聲明。按值傳遞是將參數(shù)值傳給堆棧，實參的值不隨形參的值的改變而改變；而按址傳遞是將地址傳給堆棧，實參的值隨形參的值的改變而改變，因此動態(tài)鏈接庫文件的函數(shù)或者過程的輸入?yún)?shù)通常按值傳遞，而輸出參數(shù)通常按址傳遞，如果傳遞參數(shù)為數(shù)組，則其傳遞方式為按址傳遞。

2.2 Vb程序調(diào)用.dll文件

Vb 程序要調(diào)用動態(tài)連接庫文件所含的函數(shù)或過程需先對其進行聲明，通常在Vb 工程的模塊添加聲明程序代碼來聲明動態(tài)鏈接庫文件。本文Vb 工程模塊里對動態(tài)鏈接庫文件Forlib.dll 文件的聲明程序代碼和對全局變量定義代碼為

"動態(tài)鏈接庫文件聲明

聲明程序代碼中雙引號內(nèi)的文件路徑為動態(tài)鏈接庫Forlib.dll 的存儲位置，該路徑必須詳細到根目錄，否則程序會提示找不到該動態(tài)鏈接庫文件；括號內(nèi)為Vb 調(diào)用Forlib.dll 文件的DRYCONTACT 過程傳遞的參數(shù)，其中ByVal 表示參數(shù)的傳遞按值傳遞，不可省略；ByRef 表示按址傳遞，是默認的參數(shù)傳遞方式，可省略不寫。Vb 聲明動態(tài)鏈接文件時不僅各個參數(shù)的傳遞方式要與Fortran 源程序規(guī)定的傳遞方式一致，而且其數(shù)據(jù)類型也要一致，Vb 主要變量類型與Fortran變量的對應(yīng)關(guān)系如表1 所示。

表1 Vb 與Fortran 主要變量類型的對應(yīng)關(guān)系

此外，程序中各輸入?yún)?shù)存儲變量類型也要與動態(tài)鏈接庫文件聲明中的變量類型一致。

動態(tài)鏈接庫文件聲明之后，就可以通過窗口后臺程序調(diào)用動態(tài)鏈接庫文件所包含的過程或者函數(shù)。本文Vb 程序調(diào)用動態(tài)鏈接庫文件的窗口后臺程序如下

數(shù)組的傳遞方式為按址傳遞，在調(diào)用時要把對應(yīng)數(shù)組的首個元素作為參數(shù)放在被調(diào)用函數(shù)或過程的對應(yīng)位置，如PR（1，1）、VMR（1，1）。

3 Vb 和Matlab 的混合編程

Vb 和Matlab 的混合編程主要有ActiveX 技術(shù)、可執(zhí)行文件技術(shù)、Matrix Vb 技術(shù)、動態(tài)數(shù)據(jù)交換（DDE）技術(shù)和COM 組件技術(shù)[20-21]。ActiveX 技術(shù)基于ActiveX 自動化協(xié)議，該協(xié)議允許某一個應(yīng)用軟件（控制端）通過ActiveX 自動化接口調(diào)用另一個應(yīng)用軟件（服務(wù)端）。而Vb 支持ActiveX 的控制端協(xié)議，Matlab 軟件支持ActiveX 的服務(wù)端協(xié)議。本文使用ActiveX 技術(shù)實現(xiàn)Vb 程序調(diào)用Matlab 繪圖函數(shù)進行涂層接觸分析結(jié)果的圖形化與可視化，其關(guān)鍵在于建立Vb 和Matlab 之間的ActiveX 自動化鏈接。Vb 和Matlab 之間的自動化鏈接通過在Vb 中創(chuàng)建Matlab 對象實現(xiàn)，具體方法如下

通過創(chuàng)建Vb 和Matlab 的自動化鏈接，Vb 程序就可以通過“Matlab.execute（）”等命令調(diào)用Matlab 的相關(guān)函數(shù)與命令。本文Vb 程序調(diào)用Matlab 繪制接觸壓力分布和次表層接觸應(yīng)力分布的后臺程序代碼詳見附錄2。

4 軟件運行演示驗證

軟件調(diào)試無誤后，生成CONTACT3D VERSION1.0.exe 可執(zhí)行文件。運行該執(zhí)行文件彈出如圖1 所示的軟件界面，按表2 所示的接觸參數(shù)和單元網(wǎng)格劃分數(shù)目輸入文本框。

表2 輸入的接觸參數(shù)和單元網(wǎng)格劃分數(shù)目

點擊開始計算按鈕，當軟件界面出現(xiàn)“計算完成”的提示框時，表明接觸分析結(jié)算已近完成。分別點擊接觸壓力和馮米斯應(yīng)力2 個命令按鈕，則軟件界面上顯示接觸壓力和馮米斯等效應(yīng)力的分布圖，軟件運行演示如圖3 所示。

圖3 軟件運行演示

為了進一步驗證程序計算結(jié)果的準確性，針對表2所示的接觸參數(shù)和單元網(wǎng)格劃分數(shù)目，應(yīng)用文獻[13-14]的涂層接觸分析Fortran 程序進行接觸計算，F(xiàn)ortran程序計算分析結(jié)果如圖4 所示，圖中ρ 為接觸壓力；PH代表以基體為材料的彈性半空間與球體Hertz 點接觸的最大Hertz 接觸壓力；aH代表以基體為材料的彈性半空間與球體Hertz 點接觸的最大Hertz 接觸半徑。從圖4 中可以看出，圖3 軟件界面上的結(jié)果與圖4 所示結(jié)果一致。

圖4 Fortran 程序計算分析結(jié)果

5 結(jié)語

本文基于Vb、Fortran 和Matlab 的混合編程，開發(fā)了單層涂層體系三維點接觸分析軟件，實現(xiàn)了接觸體幾何尺寸、材料特性以及接觸載荷等接觸參數(shù)輸入界面化和接觸壓力與次表層馮米斯等效應(yīng)力計算結(jié)果的界面可視化。研究表明，通過混合編程可以綜合利用Vb、Fortran 和Matlab 3 個軟件語言在界面開發(fā)、數(shù)值計算以及繪圖方面的各自優(yōu)勢，實現(xiàn)參數(shù)輸入界面化和計算結(jié)果查看界面可視化的功能，有效解決單一Fortran 代碼程序參數(shù)輸入和計算結(jié)果查看不便的問題。通過算例計算演示及與已發(fā)表論文的Fortran 程序計算結(jié)果的對比，驗證了本文開發(fā)的軟件具備預期功能，計算結(jié)果準確有效。本研究不僅為單層涂層體系接觸問題提供了一款界面友好的專用分析軟件，也為后續(xù)進一步開發(fā)多層涂層體系接觸分析軟件奠定了軟件框架和技術(shù)基礎(chǔ)。

附錄1

編譯dll 動態(tài)鏈接庫文件的Fortran 主程序文件DRYCONTACT.f90 部分程序代碼如下

附錄2

Vb 程序調(diào)用Matlab 繪制接觸壓力分布和次表層接觸應(yīng)力分布的后臺程序代碼