楊 東 生， 張 盛， 李 云 鵬＊， 陳 飆 松， 張 洪 武

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.大連理工大學(xué) 運(yùn)載工程與力學(xué)學(xué)部 工程力學(xué)系，遼寧 大連 116024)

0 引 言

工程中存在許多傳熱問題，如耐火材料的熱強(qiáng)度、極端熱環(huán)境下結(jié)構(gòu)熱防護(hù)、激光脈沖相關(guān)的高熱流及瞬態(tài)沖擊、航空航天結(jié)構(gòu)熱變形熱防護(hù)、電子電路的熱疲勞、微電子封裝以及制造業(yè)中焊接帶來的殘余應(yīng)力、工業(yè)生產(chǎn)中的鍛造成型等[1].近幾十年來，有限元法的發(fā)展為通過數(shù)值模擬解決這些工程問題提供了有效途徑[2].然而，隨著模型復(fù)雜性增加，對(duì)有限元分析程序的基本需求就是其開放性、擴(kuò)展性及大規(guī)模分析能力，而現(xiàn)有的有限元程序多基于結(jié)構(gòu)化語(yǔ)言，難以擴(kuò)展和繼續(xù)維護(hù)，因而很難滿足這樣的需求，使其進(jìn)一步發(fā)展遇到瓶頸.

面向?qū)ο蟮乃枷肟朔私Y(jié)構(gòu)化程序語(yǔ)言的不足，并為有限元程序研發(fā)提供了新的手段.面向?qū)ο蟮挠邢拊幊套钤缬?Fenves[3]提出，F(xiàn)orde等[4]給出了詳細(xì)的有限元程序的面向?qū)ο髨?zhí)行過程，文獻(xiàn)[5－7]也使用面向?qū)ο蟮木幊虂砀纳朴邢拊绦虻撵`活性.但面向?qū)ο笥邢拊幊痰难芯慷嗉杏诙x模型類，如單元、節(jié)點(diǎn)、邊界條件、材料，而很少注意和分析相關(guān)的計(jì)算任務(wù)類的定義.在模型類和分析代碼之間，沒有明顯的界限，這易使一個(gè)有限元程序變得龐大而難以維護(hù).另一方面，CAE技術(shù)是國(guó)家裝備制造業(yè)的核心競(jìng)爭(zhēng)能力之一，然而主流CAE軟件產(chǎn)品都來自國(guó)外.

面向工程與科學(xué)計(jì)算的集成化軟件平臺(tái)SiPESC[8－9]就是在這樣的背景和需求下提出來的.本文在SiPESC的基礎(chǔ)上，嘗試研發(fā)傳熱問題的面向?qū)ο笥邢拊治龀绦蛳到y(tǒng)SiPESC.THERMAL，通過插件及擴(kuò)展的管理機(jī)制，采用面向?qū)ο蟮腃＋＋程序設(shè)計(jì)語(yǔ)言、Factory及Builder等軟件設(shè)計(jì)模式，以使系統(tǒng)具有良好的開放性、擴(kuò)展性、模塊化等特點(diǎn).

本文重點(diǎn)介紹熱單元計(jì)算模塊的設(shè)計(jì)模式、單元數(shù)據(jù)類及單元計(jì)算實(shí)現(xiàn)類，以及熱荷載處理模塊的設(shè)計(jì)模式、荷載數(shù)據(jù)類及荷載處理實(shí)現(xiàn)類，以期為解決工程實(shí)際問題提供參考.

1 傳熱問題控制方程

Fourier定律是傳統(tǒng)導(dǎo)熱理論的基礎(chǔ)，其推導(dǎo)獲得Fourier傳熱控制方程[1]

其中k為傳熱系數(shù)張量，Q為熱源強(qiáng)度，T為溫度場(chǎng)，ρ為密度，c為比熱容.本文程序考慮三類熱邊界條件，即指定溫度、熱流和對(duì)流.忽略時(shí)間效應(yīng)，方程(1)退化為穩(wěn)態(tài)傳熱控制方程，即

由于Fourier理論會(huì)導(dǎo)致熱量傳播速度無限大，為了對(duì)瞬態(tài)熱傳導(dǎo)問題做出更合理的描述，需要對(duì)Fourier定律進(jìn)行修正，本文采用CV波模型，其推導(dǎo)獲得非Fourier傳熱控制方程[1]

其中τ為熱流松弛時(shí)間，其物理意義為溫度場(chǎng)的重新建立滯后于熱擾動(dòng)改變的時(shí)間.該方程為有阻尼的波動(dòng)方程，屬于雙曲型偏微分方程.

采用Galerkin離散得到穩(wěn)態(tài)傳熱、Fourier及非Fourier瞬態(tài)傳熱有限元半離散方程分別為

其中為節(jié)點(diǎn)溫度向量，對(duì)應(yīng)的整體矩陣分別由單元級(jí)別矩陣集合而成，即

單元熱傳導(dǎo)剛度、對(duì)流剛度、比熱容和熱松弛矩陣以及對(duì)流、熱流和熱源荷載矢量分別由以下積分得到:

其中N和B分別為形函數(shù)及其梯度矩陣.

對(duì)瞬態(tài)傳熱方程進(jìn)行時(shí)間離散，對(duì)Fourier和非Fourier傳熱方程分別采用θ積分[2]和Newmark直接積分法[2]進(jìn)行計(jì)算.

2 程序框架

本文程序系統(tǒng)是在SiPESC[8]及其開放式結(jié)構(gòu)有限元分析子系統(tǒng)[9]的基礎(chǔ)上完成的.SiPESC有大規(guī)模數(shù)值模型的管理能力、多物理場(chǎng)多相多尺度的集成計(jì)算能力以及系統(tǒng)開放性的特點(diǎn).

開放式結(jié)構(gòu)有限元分析子系統(tǒng)的數(shù)據(jù)存儲(chǔ)基于大規(guī)模工程數(shù)據(jù)庫(kù)系統(tǒng)，可滿足各種數(shù)據(jù)對(duì)象的設(shè)計(jì)和高效存取；數(shù)值計(jì)算部分采用模塊化設(shè)計(jì)，提供了節(jié)點(diǎn)排序、約束處理、局部坐標(biāo)轉(zhuǎn)換、單元?jiǎng)偠扔?jì)算、荷載計(jì)算、求解器等通用計(jì)算模塊.

本文程序研發(fā)在此基礎(chǔ)上展開.傳熱問題的標(biāo)量場(chǎng)特點(diǎn)，以及邊界條件的多樣性，使其具有和結(jié)構(gòu)有限元不同的特點(diǎn).因此從通用的標(biāo)量場(chǎng)問題求解角度考慮，必須設(shè)計(jì)適合標(biāo)量場(chǎng)問題求解的計(jì)算流程，如標(biāo)量場(chǎng)自由度處理、標(biāo)量場(chǎng)約束處理等，這些流程不但適合傳熱問題的分析，也適合其他標(biāo)量場(chǎng)如滲流場(chǎng)等問題的處理.

考慮到傳熱問題的本構(gòu)、荷載的特殊性，本文增加了傳熱問題特有的計(jì)算流程，如單元矩陣的計(jì)算包括單元熱傳導(dǎo)陣、單元比熱容陣，以及后處理流程單元溫度梯度、單元熱流及節(jié)點(diǎn)熱流計(jì)算，程序框架如圖1所示.

圖1 分析系統(tǒng)程序框架Fig.1 Program framework of the analysis system

各個(gè)流程之間采用數(shù)據(jù)定義標(biāo)準(zhǔn)接口，耦合性大大降低，各流程分別調(diào)用不同的插件完成各自的功能，對(duì)不同問題可以通過不同的流程動(dòng)態(tài)組裝得到不同的分析類型.穩(wěn)態(tài)熱分析的流程為節(jié)點(diǎn)排序、標(biāo)量場(chǎng)自由度處理、標(biāo)量場(chǎng)約束處理、單元熱傳導(dǎo)陣計(jì)算、荷載計(jì)算及求解.由于定義了標(biāo)準(zhǔn)的數(shù)據(jù)接口，不同的分析類型可以共用相同的流程，體現(xiàn)了模塊化的靈活性.

3 熱單元和熱荷載計(jì)算模塊

熱單元和熱荷載計(jì)算是本文程序系統(tǒng)的重要計(jì)算流程，本文將介紹其實(shí)現(xiàn)過程.為了提高代碼的可重用性、提高軟件質(zhì)量、縮短研發(fā)時(shí)間，在熱單元和熱荷載計(jì)算模塊采用了軟件設(shè)計(jì)模式中的Factory模式和Builder模式[10]，使本文程序系統(tǒng)具有良好的開放性、擴(kuò)展性、模塊化等特點(diǎn).

3.1 Factory模式

在有限元計(jì)算時(shí)，對(duì)不同單元需生成不同的計(jì)算模塊，傳統(tǒng)方法往往用new生成.以荷載計(jì)算類為例，假設(shè)要滿足節(jié)點(diǎn)熱荷載及單元熱流荷載的計(jì)算，傳統(tǒng)的面向?qū)ο蠓椒ㄊ紫榷x一個(gè)荷載計(jì)算基類，由此類分別派生出節(jié)點(diǎn)熱荷載計(jì)算類及單元熱流荷載計(jì)算類，其類圖如圖2(a)所示.則對(duì)節(jié)點(diǎn)熱荷載的調(diào)用方式如下:

LoadCalculator calculator＝

new NodeHeatLoadCalculator；

loadComponent＝calculator.start(load)；

對(duì)單元熱流荷載的調(diào)用方式如下:

LoadCalculator calculator＝

new EleHeatFluxLoadCalculator；

loadComponent＝calculator.start(load)；

可以看出，對(duì)不同種類的荷載，調(diào)用方式不同，這樣添加新的荷載類型時(shí)需要修改已有的調(diào)用代碼，降低了程序的可維護(hù)性和可擴(kuò)展性.

本文采用軟件設(shè)計(jì)模式中的Factory模式，對(duì)象工廠首先生成父工廠類，定義通過單元類型變量創(chuàng)建對(duì)象的公共接口，而其子類則負(fù)責(zé)生成具體的類實(shí)例.這樣可以不用關(guān)心具體對(duì)象的實(shí)現(xiàn)，并可通過工廠提供的優(yōu)先級(jí)方法動(dòng)態(tài)替換單元計(jì)算模塊，提高了系統(tǒng)的靈活性.

圖2 荷載計(jì)算類圖Fig.2 Class diagram for load calculation

在具體實(shí)施過程中，對(duì)不同的單元類型(或荷載、積分、形函數(shù))定義不同的type，在計(jì)算時(shí)由此type找到對(duì)應(yīng)的Factory，然后由Factory創(chuàng)建出對(duì)應(yīng)的計(jì)算模塊.

仍以荷載計(jì)算類為例，本文方法首先定義荷載計(jì)算基類與荷載計(jì)算工廠基類，由此派生出不同種類的荷載計(jì)算類及荷載計(jì)算工廠類.其類圖如圖2(b)所示.對(duì)應(yīng)的調(diào)用方式如下:

type＝load.getType()；

MLoadCalculatorFactory factory＝

factoryManager.getFactory(type)；

MLoadCalculator calculator＝factory.create()；

loadComponent＝calculator.start(load)；

即首先由荷載類型找到對(duì)應(yīng)的荷載計(jì)算工廠，然后創(chuàng)建出對(duì)應(yīng)的計(jì)算類，統(tǒng)一完成計(jì)算.由于調(diào)用方式統(tǒng)一，添加新荷載不需要修改已有的調(diào)用代碼，方便維護(hù)和擴(kuò)展.

3.2 Builder模式

在單元矩陣的計(jì)算過程中需要多次數(shù)值積分，如式(6)所示.數(shù)值積分算法較穩(wěn)定，但是單元類型和材料模型的多樣性使穩(wěn)定的數(shù)值積分算法不能重用.本文采用軟件設(shè)計(jì)中的Builder模式，采用封裝機(jī)制來隔離出復(fù)雜對(duì)象的各個(gè)部分的變化，保持系統(tǒng)中穩(wěn)定構(gòu)建算法不隨其子對(duì)象的改變而改變，使同樣的構(gòu)造過程可創(chuàng)建不同的表示.以單元?jiǎng)偠汝?單元熱傳導(dǎo)陣)積分計(jì)算為例，在單元計(jì)算前，由MElementStiffBuilder對(duì)象將所需對(duì)象生成，傳遞給單元計(jì)算類 MElementStiff－Calculator，類圖如圖3所示.

圖3 Builder模式的單元?jiǎng)偠汝囉?jì)算Fig.3 Builder pattern for element stiffness matrix calculation

這樣單元計(jì)算對(duì)象就只負(fù)責(zé)計(jì)算任務(wù)，而不必判斷單元的性質(zhì).對(duì)不同類型的單元積分，只需定義新的構(gòu)造器，例如三維六面體單元熱傳導(dǎo)陣積分設(shè)置為

builder.setConstiType("heatbrickeleconstidparser")；

builder.setShpGradType(HexaBrickShpGradMatrix)；

builder.setIntegType(Gau3D2IntegForm)；

三維四面體單元熱傳導(dǎo)陣積分設(shè)置為

builder.setConstiType("heatbrickeleconstidparser")；

builder.setShpGradType(TetraBrickShpGradMatrix)；

builder.setIntegType(TetraIntegForm)；

對(duì)于不同的單元，可以重用計(jì)算模塊，體現(xiàn)了軟件設(shè)計(jì)的可重用性.

3.3 熱單元計(jì)算模塊

工程中傳熱問題的多樣性會(huì)導(dǎo)致有限元分析單元類型多種多樣，一些通用有限元軟件甚至有上百種單元.這些單元通常還不能滿足用戶的需求，隨著理論和算法的發(fā)展，會(huì)出現(xiàn)更多的單元類型.所以對(duì)于一個(gè)開放式有限元計(jì)算平臺(tái)，一個(gè)重要功能是可以提供開放的接口，方便添加新的單元.

為了提高系統(tǒng)的靈活性，本文數(shù)據(jù)存儲(chǔ)采用工程數(shù)據(jù)庫(kù)系統(tǒng)，可滿足各種數(shù)據(jù)對(duì)象的設(shè)計(jì)和高效存取.計(jì)算部分采用模塊化設(shè)計(jì)，各個(gè)流程間采用數(shù)據(jù)定義標(biāo)準(zhǔn)接口，通過工程數(shù)據(jù)庫(kù)系統(tǒng)完成各自功能.本文單元計(jì)算模塊采用軟件設(shè)計(jì)模式中的Factory和Builder模式，建立通用單元計(jì)算模塊.

(1)單元數(shù)據(jù)類

對(duì)于不同類型的單元，首先需要定義單元數(shù)據(jù)類.單元數(shù)據(jù)類由工程數(shù)據(jù)庫(kù)系統(tǒng)的數(shù)據(jù)基類MDataObject派生而來，MDataObject類只定義數(shù)據(jù)對(duì)象ID的存取方法.為了保證各單元數(shù)據(jù)使用過程中的通用性，定義單元數(shù)據(jù)基類MElementData，該基類定義的主要接口為獲取單元類型、存取單元指定序號(hào)的節(jié)點(diǎn)ID、獲取單元節(jié)點(diǎn)個(gè)數(shù)、存取單元性質(zhì)ID、獲取節(jié)點(diǎn)剛度自由度貢獻(xiàn)描述數(shù)組.具體某種單元的數(shù)據(jù)類可由上述基類派生，以標(biāo)量場(chǎng)單元數(shù)據(jù)為例，單元基類接口定義以及數(shù)據(jù)類之間的派生關(guān)系如圖4所示.

圖4 數(shù)據(jù)基類、單元數(shù)據(jù)基類及單元數(shù)據(jù)類類圖Fig.4 Class diagram for data base class，element data base class and element data classes

對(duì)某種單元，單元類型、單元節(jié)點(diǎn)個(gè)數(shù)、自由度貢獻(xiàn)描述數(shù)組已確定，只需從數(shù)據(jù)庫(kù)存取指定位置節(jié)點(diǎn)ID和單元性質(zhì)ID，這大大簡(jiǎn)化了單元數(shù)據(jù)類的設(shè)計(jì).每個(gè)單元數(shù)據(jù)類均有唯一的類型type作為標(biāo)識(shí)，在各計(jì)算流程中，若計(jì)算過程根據(jù)各單元不同，則通過該type和工廠模式可找到對(duì)應(yīng)的計(jì)算模塊.

(2)單元計(jì)算實(shí)現(xiàn)類

在SiPESC結(jié)構(gòu)有限元[8]系統(tǒng)中，已定義了應(yīng)變陣、本構(gòu)解析、材料解析等基類和接口，本文的形函數(shù)梯度陣、熱傳導(dǎo)本構(gòu)解析、傳熱材料解析均由這些基類派生而來.

SiPESC的插件和擴(kuò)展機(jī)制定義了擴(kuò)展基類MExtension，所有擴(kuò)展均由此派生.為了更好地管理和分配不同層次的計(jì)算任務(wù)，在擴(kuò)展基類的基礎(chǔ)上，定義了三種基類，即任務(wù)基類MTask、任務(wù)執(zhí)行基類MTaskExecutor和任務(wù)管理基類MTaskManager.其派生關(guān)系類圖如圖5所示.

圖5 擴(kuò)展基類、任務(wù)基類、任務(wù)執(zhí)行基類、任務(wù)管理基類類圖Fig.5 Class diagram for extension base class，task base class，task executor base class and task manager base class

其中MTask類定義的虛方法initialize完成數(shù)據(jù)庫(kù)的初始化.MTaskExecutor類定義的虛方法start完成單一任務(wù)(比如單一單元、單一荷載)的計(jì)算，其輸入輸出參數(shù)均為以MDataObject為基類的數(shù)據(jù)類.MTaskManager類定義的虛方法start打開數(shù)據(jù)庫(kù)，把任務(wù)分配給 MTaskExecutor，得到計(jì)算結(jié)果，并存入數(shù)據(jù)庫(kù).

單元計(jì)算模塊采用Factory模式和Builder模式結(jié)合，其中Factory基類MExtensionFactory由擴(kuò)展基類MExtension派生而來，Builder基類MElementBuilder由MTask派生而來，圖6給出了各個(gè)實(shí)現(xiàn)類之間的關(guān)系.

其中Factory基類MExtensionFactory定義的成員函數(shù)有g(shù)etType和createExtension，其中前者返回單元類型type，供擴(kuò)展工廠管理器查找，后者調(diào)用MHeatEleStiffBuilder的3個(gè)成員函數(shù)，完成對(duì)本構(gòu)、形函數(shù)及其梯度、積分類型的設(shè)置之后，返回 MHeatEleStiffBuilder對(duì)象.而MHeatEleStiffBuilder分別由Factory模式創(chuàng)建出對(duì)應(yīng)的計(jì)算模塊，通過調(diào)用 MHeatEleStiff－Calculator的3個(gè)成員函數(shù)，完成計(jì)算模塊的設(shè)置，通過成員函數(shù)build返回設(shè)置好的MHeatEle－StiffCalculator.最 后 由 MHeatEle－StiffCalculator繼承基類MTaskExecutor的start方法完成對(duì)式(6a)的積分計(jì)算.

圖6 傳熱單元?jiǎng)偠汝囉?jì)算類圖Fig.6 Class diagram for thermal element stiffness matrix calculation

對(duì)傳熱桿單元，其單元?jiǎng)偠汝囉上率浇馕鼋o出:

故不需要采用上述積分方法計(jì)算，但為了各個(gè)單元的統(tǒng)一，仍采用同樣的流程，只是不用設(shè)置本構(gòu)、形函數(shù)及其梯度、積分類型等信息，如圖6所示.

值得注意的是，所有單元模塊的計(jì)算對(duì)象都由Factory模式生成，這樣各計(jì)算對(duì)象都可以靈活地替換和組裝.這樣對(duì)于不同的單元，只需完成對(duì)應(yīng)的本構(gòu)、形函數(shù)及其梯度、積分類型等模塊的定義，即可重用已有的單元?jiǎng)偠汝囉?jì)算模塊.

3.4 熱荷載處理模塊

在傳熱問題的有限元分析中，荷載形式多種多樣.根據(jù)荷載類型不同，可以分為熱源荷載、熱流荷載、對(duì)流荷載等.根據(jù)施加荷載幾何形式不同，可分為體荷載、面荷載、點(diǎn)荷載.根據(jù)施加有限元對(duì)象不同，可以分為節(jié)點(diǎn)荷載和單元表面荷載.

在不同荷載的處理上采用Factory模式，即每種荷載均含有一個(gè)荷載類型type，根據(jù)該type找到對(duì)應(yīng)的荷載計(jì)算器進(jìn)行相應(yīng)的處理.

對(duì)于和單元無關(guān)的荷載，如節(jié)點(diǎn)熱源，處理較簡(jiǎn)單.對(duì)于和單元相關(guān)的荷載，多樣的單元類型和荷載形式使開發(fā)者的工作變得很煩瑣.比如開發(fā)者若想增加一種單元，則在完成單元計(jì)算的同時(shí)，需要處理此種單元表面施加的不同荷載形式.若新增加一種和單元相關(guān)的荷載，則需要對(duì)所有支持此種荷載的單元開發(fā)出對(duì)應(yīng)的計(jì)算模塊.這樣使開發(fā)的耦合度上升，且不利于代碼的重用性.本文程序系統(tǒng)采用荷載面單元的處理方式，有效地降低了兩者之間的耦合度，使開發(fā)的工作量降低.

(1)荷載數(shù)據(jù)類

以表面荷載(熱流或?qū)α?為例，不同幾何形式表示的表面荷載可分為面、線(需定義厚度)、點(diǎn)(需定義面積)3種形式.其中每種幾何形式可以有不同的類型，比如面形式也可以是三節(jié)點(diǎn)三角形或六節(jié)點(diǎn)三角形、八節(jié)點(diǎn)四邊形等高階形式.

對(duì)不同的幾何形式，需定義幾何表面單元數(shù)據(jù)類作為一種新的單元類型，但是不需要定義自由度信息，為了方便統(tǒng)一處理，定義幾何表面單元數(shù)據(jù)基類MGeoFaceData.對(duì)應(yīng)的荷載類只需定義荷載類型、荷載值和該荷載施加幾何對(duì)象的ID.

在表面荷載的處理上，采用Factory的設(shè)計(jì)模式，即根據(jù)幾何表面單元的type找到對(duì)應(yīng)的幾何表面單元荷載計(jì)算器，完成計(jì)算并返回各個(gè)節(jié)點(diǎn)的荷載值.

(2)荷載處理實(shí)現(xiàn)類

對(duì)于面荷載，通常由各個(gè)節(jié)點(diǎn)描述一個(gè)荷載密度曲面，即定義各個(gè)節(jié)點(diǎn)的荷載密度值，曲面上的其他點(diǎn)由一定的插值方式得到.此種形式荷載由數(shù)值積分轉(zhuǎn)化為節(jié)點(diǎn)荷載.

對(duì)于點(diǎn)形式的幾何表面單元，其節(jié)點(diǎn)荷載的轉(zhuǎn)化可由解析形式完成，即

其中f為表面荷載密度，A為荷載施加面積，F(xiàn)為節(jié)點(diǎn)荷載向量.二節(jié)點(diǎn)線幾何表面單元的積分也有解析解，其荷載項(xiàng)和對(duì)流剛度項(xiàng)分別為

其中L為線長(zhǎng)度，t為厚度，f1和f2分別為兩個(gè)節(jié)點(diǎn)上荷載密度大小，兩點(diǎn)之間的荷載密度由這兩個(gè)節(jié)點(diǎn)的值插值得到.

在程序?qū)崿F(xiàn)上采用Factory模式，類圖如圖7所示.由節(jié)點(diǎn)荷載、對(duì)流荷載、熱流荷載的type找到對(duì)應(yīng)的Factory，由Factory創(chuàng)建出對(duì)應(yīng)的計(jì)算模塊.其中各對(duì)應(yīng)模塊均有共同的基類MLoad－Calculator，此基類繼承于 MTaskExecutor任務(wù)執(zhí)行基類.

圖7 節(jié)點(diǎn)荷載、對(duì)流荷載和熱流荷載類圖Fig.7 Class diagram for node load，convection load and heat flux load

其中節(jié)點(diǎn)荷載計(jì)算與單元無關(guān)，直接由MScalarNodeLoadCalculator完成，而對(duì)流荷載和熱流荷載需要對(duì)單元進(jìn)行面積分，所以分別調(diào)用對(duì)應(yīng)的積分計(jì)算模塊完成.其中式(6e)～(6g)的荷載面積分由MScalarFaceLoadCalculator模塊完成，式(6b)中的對(duì)流剛度積分由MScalarFaceStiffCalculator模塊完成.這兩個(gè)模塊的積分與單元?jiǎng)偠饶K的積分類似，采用Factory模式與Builder模式相結(jié)合的設(shè)計(jì)模式，如圖8所示.

本文考慮了三節(jié)點(diǎn)三角形、四節(jié)點(diǎn)四邊形，及二節(jié)點(diǎn)線、點(diǎn)等不同形式的表面單元，對(duì)于需要數(shù)值積分的面單元，由Builder模式、Factory模式中給定的type參數(shù)創(chuàng)建出插值形函數(shù)對(duì)象、形函數(shù)梯度對(duì)象、積分格式對(duì)象，傳遞給Calculator，統(tǒng)一完成計(jì)算.對(duì)于有解析解的點(diǎn)、線等幾何形式的表面單元，仍然采用統(tǒng)一的流程，由對(duì)應(yīng)的Calculator直接計(jì)算完成.

圖8 熱荷載積分和熱對(duì)流剛度積分類圖Fig.8 Class diagram for heat load integration and heat convection stiffness integration

4 數(shù)值算例

本節(jié)給出了一維、二維、三維常用單元在指定溫度、對(duì)流、熱流等不同邊界條件下的算例驗(yàn)證.

4.1 一維散熱片問題

矩形鋁制散熱片[11]，表面溫度為100℃，環(huán)境溫度20℃，鋁熱導(dǎo)率為168W／(m·K)，對(duì)流傳熱系數(shù)為30W／(m2·K)，散熱片長(zhǎng)80mm，寬5mm，厚1mm.采用一維單元建模，劃分5個(gè)節(jié)點(diǎn)，4個(gè)熱傳導(dǎo)單元，4個(gè)對(duì)流邊界單元，1個(gè)對(duì)流邊界節(jié)點(diǎn)(節(jié)點(diǎn)5)，如圖9所示.

圖9 一維散熱片模型與單元離散Fig.9 Model and mesh for one－dimension cooling plate

計(jì)算得到的結(jié)果與文獻(xiàn)[11]中的結(jié)果比較如表1所示，從表中可以看出與文獻(xiàn)結(jié)果吻合較好.

4.2 二維散熱片問題

考慮電子設(shè)備中的鋁制散熱片[11]，熱導(dǎo)率為170W／(m·K)，基座熱流密度為1 000W／m2，環(huán)境溫度為20℃，對(duì)流傳熱系數(shù)為40 W／(m2·K).

由于對(duì)稱性，取一半進(jìn)行建模，劃分為四邊形網(wǎng)格，共有3 093個(gè)節(jié)點(diǎn)，2 743個(gè)單元.計(jì)算得到的溫度分布云圖如圖10所示.

從圖中可以看出，結(jié)構(gòu)溫度分布較合理，溫度峰值為48.463℃，與ANSYS軟件計(jì)算值的前5位完全相同.

表1 一維散熱片溫度分布結(jié)果比較Tab.1 Temperature distribution results comparison of one－dimension cooling plate

圖10 二維散熱片溫度分布云圖Fig.10 Temperature contour plot for two－dimension cooling plate

4.3 換熱管道傳熱

為了提高換熱管道的傳熱特性，在內(nèi)壁形成一系列軸向散熱片，其中材料的熱導(dǎo)率為400W／(m·K)，內(nèi)部流體溫度80℃，內(nèi)壁對(duì)流傳熱系數(shù)為150W／(m2·K)，外部環(huán)境溫度15℃，外壁對(duì)流傳熱系數(shù)30W／(m2·K).

劃分為三角形網(wǎng)格，共有3 632個(gè)節(jié)點(diǎn)，6 484個(gè)單元，計(jì)算得到的溫度分布云圖如圖11所示.

圖11 換熱管道溫度分布云圖Fig.11 Temperature contour plot for heat exchange tube

從圖中可以看出，管道溫度分布從72.806～73.500℃，與ANSYS軟件計(jì)算值的前5位完全相同.

4.4 齒輪熱分析

齒輪運(yùn)轉(zhuǎn)過程中，由于齒間摩擦生熱會(huì)對(duì)齒輪的性能產(chǎn)生影響，本算例模擬齒輪運(yùn)轉(zhuǎn)過程中的溫度分布情況，其中材料熱導(dǎo)率為230 W／(m·K)，內(nèi)壁固定溫度0℃，各齒緣與另一齒輪接觸處給定熱流密度是100W／m2來模擬運(yùn)轉(zhuǎn)過程中的摩擦生熱情況.

劃分六面體網(wǎng)格，其中有55 910個(gè)節(jié)點(diǎn)，46 200個(gè)單元，計(jì)算得到的溫度分布云圖如圖12所示.

圖12 齒輪運(yùn)轉(zhuǎn)過程中溫度分布云圖Fig.12 Temperature contour plot for working gear

從圖中可以看出，溫度峰值為16.510℃，與ANSYS軟件計(jì)算值的前5位完全相同.

5 結(jié) 語(yǔ)

本文基于SiPESC集成平臺(tái)，設(shè)計(jì)了傳熱問題有限元分析的通用程序系統(tǒng)，從單元計(jì)算模塊和荷載計(jì)算模塊兩方面介紹了程序系統(tǒng)的設(shè)計(jì)模式.本文程序基于SiPESC的插件和擴(kuò)展機(jī)制，采用C＋＋面向?qū)ο笤O(shè)計(jì)語(yǔ)言及軟件設(shè)計(jì)模式中的Factory模式和Builder模式，程序特點(diǎn)如下:

(1)代碼重用.提取相同算法作為通用模塊，比如對(duì)不同單元矩陣的數(shù)值積分算法進(jìn)行提取，以實(shí)現(xiàn)代碼重用.

(2)開放性接口.程序所有接口全部開放，便于添加新的單元、荷載和功能.

(3)良好的可維護(hù)性.由于采用插件和擴(kuò)展機(jī)制，以及Factory設(shè)計(jì)模式，添加新的單元或荷載不需對(duì)原有代碼進(jìn)行修改，程序原有功能也可由Factory的優(yōu)先級(jí)管理實(shí)現(xiàn)動(dòng)態(tài)替換.

本文程序系統(tǒng)有良好的開放性和擴(kuò)展性，可在此基礎(chǔ)上進(jìn)一步展開瞬態(tài)傳熱、非線性傳熱以及熱力耦合等問題的程序研發(fā)，也可以用于其他標(biāo)量場(chǎng)相關(guān)問題的分析和研究.

[1]Tamma K K， Namburu R R.Computational approaches with applications to non－classical and classical thermomechanical problems [J].Applied Mechanics Reviews，1997，50(9):514－551.

[2]Zienkiewicz O C，Taylor R L，Zhu J Z.The Finite Element Method:Its Basis and Fundamentals[M].6th ed.Oxford:Elsevier Butterworth－Heinemann，2005.

[3]Fenves G L. Object－oriented programming for engineering software development[J].Engineering with Computers，1990，6(1):1－15.

[4]Forde B W R，F(xiàn)oschi R O，Stiemer S F.Objectoriented finite element analysis [J].Computers and Structures，1990，34(3):355－374.

[5]Archer G C，F(xiàn)enves G，Thewalt C.A new objectoriented finite element analysis program architecture[J].Computers and Structures，1999，70(1):63－75.

[6]張 向，許晶月，沈啟彧，等.面向?qū)ο蟮挠邢拊绦蛟O(shè)計(jì)[J].計(jì)算力學(xué)學(xué)報(bào)，1999，16(2):216－226.ZHANG Xiang，XU Jing－yue，SHEN Qi－yu，etal.Object－oriented finite element programming [J].Chinese Journal of Computational Mechanics，1999，16(2):216－226.(in Chinese)

[7]馬永其，馮 偉.面向?qū)ο笥邢拊治龀绦蛟O(shè)計(jì)及其VC＋＋實(shí)現(xiàn)[J].應(yīng)用數(shù)學(xué)和力學(xué)，2002，23(12):1283－1288.MA Yong－qi，F(xiàn)ENG Wei.Object－oriented finite element analysis and programming in VC＋＋ [J].Applied Mathematics and Mechanics，2002，23(12):1283－1288.(in Chinese)

[8]張洪武，陳飆松，李云鵬，等.面向集成化CAE軟件開發(fā)的SiPESC研發(fā)工作進(jìn)展[J].計(jì)算機(jī)輔助工程，2011，20(2):39－49.ZHANG Hong－wu，CHEN Biao－song，LI Yunpeng，etal. Advancement of design and implementation of SiPESC for development of integrated CAE software systems [J].Computer Aided Engineering，2011，20 (2):39－49. (in Chinese)

[9]張 盛，楊東生，尹 進(jìn)，等.SiPESC.FEMS的單元計(jì)算模塊設(shè)計(jì)模式[J].計(jì)算機(jī)輔助工程，2011，20(3):49－54.ZHANG Sheng，YANG Dong－sheng，YIN Jin，etal. Design pattern of element computation module of SiPESC.FEMS [J].Computer Aided Engineering，2011，20(3):49－54.(in Chinese)

[10]Gamma E，Helm R，Johnson R，etal.Design Patterns:Elements of Reusable Object－oriented Software[M].Boston:Addison Wesley，1995:11－45.

[11]Moaveni S.有限元分析——ANSYS理論與應(yīng)用[M].3版.王 崧，劉麗娟，董春敏，等，譯.北京:電子工業(yè)出版社，2008.Moaveni S.Finite Element Analysis— Theory and Application with ANSYS [M].3rd ed.WANG Song，LIU Li－juan，DONG Chun－min，etal，tran.Beijing:Publishing House of Electronics Industry，2008.(in Chinese)