三維對流擴散問題有限元計算

摘要對流擴散問題的數(shù)值求解是傳熱、傳質學的重要研究課題。研究對流擴散問題的數(shù)值解，尋找一種快速、穩(wěn)定而又實用的數(shù)值算法，具有重要的理論和實際意義。目前，用來求解這…類問題的數(shù)值解法有很多，研究的角度也有所不同。在各種數(shù)值方法的求解中，其程序開發(fā)的具體實現(xiàn)一般都采用結構化語言，用面向對象的程序設計語言來開發(fā)的相對較少。本文研究的課題是面向對象的三維對流擴散問題有限元分析程序的設計與實現(xiàn)，主要研究Peclet數(shù)較小時，三維定常的對流擴散方程的有限元解法，用面向對象的C++程序設計語言開發(fā)出有應用價值的通用程序，對相關的對流擴散問題進行數(shù)值模擬。本文的主要內容總結如下：.闡述了本文的研究背景，綜述了面向對象程序設計的基本原理。.采用四結點四面體單元和八結點六面體等參數(shù)單元進行三維定常對流擴散方程的有限元分析。.用面向對象的程序設計語言C++來具體實現(xiàn)三維定常對流擴散問題的有限元程序的開發(fā)，融入用C++編制的大型通用有限元科學和工程計算軟件。.給出了四個算例，并對分析結果進行初步討論和評價，從而驗證本文方法的有效性和實用性。關鍵詞：對流擴散方程，有限元，面向對象方法，程序設計ABSTRACTNumericalsolutionofconvection-diffusionproblemisanimportantsubjectinheatandmasstransfer.Tothisproblem,itprovidesimportanttheoreticandpracticalsignificancetoseekafast,stableandpracticalnumericalmethod.Atpresent,therearemanymethodstosolvethisclassofproblems,withitsdiflferentpoints.Amongthesemethods,designandimplementationofprogramsweredevelopedbystructurelanguagesmostly,andlessbyobject-orientedprogramcomparatively.Thesubjectofthispaperistodevelopanobject-orientedfiniteelementanalysisprogramofthree-dimensionalconvection-diffusionequation.Wemainlystudythefiniteelementmethodforthree-dimensionalsteady-stateconvection-diffusionequationinthecaseofthevalueofPecletissmall,andthendevelopageneralprogramwithobject-orientedC++languageandgivethenumericalsimulationtothecorrespondingconvection-diffusionproblems.Themainofthispaperisoutlinedasfollows:.Thebackgroundofthispaperisdescribed,andthebasictheoryofobject-orientedprogramissummarized..Thefiniteelementanalysisofthree-dimensionalsteady-stateconvection-diffusionequationisgivenbyusing4-nodetetrahedralelementsand8-nodehexahedralelements..Thefiniteelementprogramofthree-dimensionalsteady-stateconvection-diffusionproblemsisdevelopedbyusingtheobject-orientedC++language,whichcanbeappliedtoalargegeneralfiniteelementscienceandengineeringcomputationsoftware..Theprimarydiscussionandtheremarksabouttheresultsoffournumericalexamplesaremade,andthenthevalidityandutilityofthemethodproposedinthispaperisverified.Keywords:Convection-diffusionequation,Finiteelementmethod,Object-orientedmethod,Programming1緒論對流擴散方程是一類基本的運動方程，它可以用來描述河流污染、大氣污染、核廢物污染中污染物質的分布，流體的流動和流體中熱的傳導等眾多物理現(xiàn)象。對于同時伴有物質輸運和分子擴散的物理過程以及粘性流體的流動，其數(shù)學模型通常為對流擴散方程或含有此類方程的偏微分方程組的定解問題。對流擴散問題數(shù)值計算方法的研究具有重要的理論和實際意義，可用于環(huán)境科學、能源開發(fā)、流體力學和電子科學等許多領域。對流擴散問題的有效數(shù)值解法一直是計算數(shù)學中重要的研究內容，求解對流擴散方程的數(shù)值方法主要是有限差分法(FDM)、有限元法(FEM)、有限體積法(FVM)、有限解析法(FAM)、邊界元法(BEM)、譜方法(SM)等多種方法。但是對于對流占優(yōu)問題，用通常的差分法或有限元法進行求解將出現(xiàn)數(shù)值震蕩。為了克服數(shù)值震蕩，80年代，J.Douglas,Jr.和T.F.Russell等提出特征修正技術求解對流擴散占優(yōu)的對流擴散問題，與其它方法相結合，提出了特征有限元方法、特征有限差分方法、特征混合元方法；T.J.Hughes和A.Brooks提出過一種沿流線方向附加人工黏性的間斷有限元法，稱為流線擴散方法(SDM)。有限差分法、有限元法、有限體積法是工程應用中的主要方法。本文所研究是低Peclet數(shù)下，粘性不可壓縮流體的流動中，三維定常的對流擴散問題的有限元解法。驗證Peclet數(shù)較小時，三維定常的對流擴散方程的有限元解法的有效性，用面向對象的C++程序設計語言開發(fā)出有應用價值的通用程序，對相關的對流擴散問題進行數(shù)值模擬。對于對流占優(yōu)對流擴散問題的求解，采用本文方法，要得到穩(wěn)定解，則要通過加密有限元網格來實現(xiàn)。對流擴散方程的背景及其特點對流擴散問題，又稱彌散(Dispersion)問題。它主要是研究流體中由流體質點所攜帶的某種物理量，如溫度或溶解于流體中的物質的濃度在流動過程中的變化規(guī)律。這些變化一般包括對流、擴散以及由于某種物理化學的原因引起的物理量自身衰減或增長的過程。對流擴散方程的背景四假定u=u(x,y,z,t)是流體中單位體積所攜帶的物理量。如物理量是流體介質的質量，則u表示密度；如物理量是某種污染物的含量，則u是表示污染濃度；如物理量是熱能，則u是熱量或溫度。為研究物理量u的變化規(guī)律，在流場中任取一個有限區(qū)域丫，其邊界為SoV中流體所攜帶的物理量u的變化，由下列三方面物理過程組成。(1)對流過程由于流體處于宏觀的流動狀態(tài)，丫中U的總量因對流而發(fā)生變化，可以認為它由兩方面組成：一項是在u中V隨著時間的增長而產生的變化；另一項是區(qū)域V中由于流體流動，位置發(fā)生變化所引起流體中u的變化(圖1.1)o也就是V中u的積分量的隨體導數(shù)V(r+Ar) SVnV(r)圖1.1對流過程Figure1.1Processofconvectiondtdu

dtdt其中Vn=v-n是流體速度在S面上的法向分量。利用Green-Gauss公式Nmv/iJS=\div(uv)dV應有dt(2)擴散過程擴散過程包含有分子布朗運動形成的分子擴散以及流體湍流運動形成的湍流擴散。這種擴散使得物理量u在流場中由高值向低值方向移動。若擴散的速率用q表示，它是單位時間內通過單位面積的物理量，F(xiàn)ick定律給出了物理量在流場中的擴散速率公式q="Kgradu其中張量K=Km+K'm的各個分量值取決于含有物理量u的流體的狀態(tài)與性質。-一般地說，企圖在理論在確定K值是非常困難的，通常是通過實驗測定。如果擴散是各向同性的，則可寫成q="kgradu其中k是擴散系數(shù)，它可以是M或其它物理量的函數(shù)，也可以是常數(shù)。由于擴散作用，使區(qū)域V中U的增加量為，K是分子擴散系數(shù)張量，K,是湍流擴散系數(shù)張量。K-N?qdS=N?KgradudS=Jdiv{Kgradu)dV(3)源和匯流場中物理量u的自身增長與衰減，一般是通過分布在流場中的源和匯來描述。這種源和匯是場的分布函數(shù)，記為。。。>0表示源，它將使u增長；。V0表示匯，它將使u減小。m增長或減小的速率由源或匯的強度，即。絕對值的大小來反映。在區(qū)域V中，由于源匯的作用使M的增加量為通常我們將Q分為源增長項和匯衰減項。線性衰減類型的表達式為其中7>0是源分布函數(shù)，PN0是衰減函數(shù)。另一種常用的衰減類型的時間的指數(shù)規(guī)律Q=f-uoe7其中wo是初始時刻物理量〃的濃度，t是時間變量，y是衰減常數(shù)。1.1.2對流擴散方程的導出?則根據(jù)守恒原理，單位時間內丫中M的變化應滿足：［由于對流M的增加量］=［由于擴散進入V中的M的增加量］+［由于源匯作用U的增加量］也就是由于V是任意的，上面的守恒方程可寫成cu一+div(uv)=div(Kgradu)+Q (Lia)a .或寫成加.74-vgradu=div(Kgradu)+Q (1.1b)立這個方程通常稱為對流擴散方程。方程(1.1)的邊界條件分為兩種類型：(1)在邊界口上給出本質邊界條件(2)在邊界n上給出自然邊界條件,3CU(Jn=Sx,rii=qnJl初+div(uv)|JdV=Jdiv(Kgradu)dV+jQdV總邊界「=BU「2。上述邊界條件也可統(tǒng)一寫成：在邊界「上SuUokij n：+a,M+az=0加式中a。、ai、a2是已知函數(shù)。初始條件為：在t=0時，給出〃的分布H(X,0)=H(X)方程(1.1)式中的速度V一般有兩種情況：一種是U對于流場中的密度P或速度V沒有直接影響，或者影響很小可以忽略，此時V可以獨立通過流體力學的主流方程(即連續(xù)性方程和運動方程)求出。另外--種是U和流場中的密度P或速度V有直接的關系，必須耦合求解。本文進行的有限元分析，考慮前面一種情形，并認為已經給出了速度場Vo另外由于本文是研究三維定常的對流擴散方程的有限元解法，則流體中所攜帶的物理量M與時間t無關，此時，在方程(1.1)中也=0

a對流擴散方程是流體力學中的典型方程，連續(xù)方程、動量方程、能量方程都可以寫成對流擴散方程的形式。1.1.3對流擴散方程的特點啪⑸對流擴散方程(Lia)和(1.1b)右端第一項為擴散項，左端第二項則是對流項。由于其方程本身的特點，給建立準確有效的數(shù)值求解方法帶來一定的困難。對流和擴散給流體中由流體攜帶的某種物理量的變化過程，可以通過一個無量綱的特征參數(shù)(Peclet數(shù))來描述，Peclet數(shù)Pe的定義為：Pe="

D這里v是來流速度，L是特征長度，D是物質的擴散系數(shù)。如果Pe數(shù)較小，即對流效應相對較弱，這類問題中，擴散占主導地位，方程(1.1a)和(1.1b)是橢圓型或拋物線型；如果Pe數(shù)較大，即溶質分子的擴散相對于流體速度而言是緩慢的，這類問題中，對流占優(yōu)，方程(LIa)和(1.1b)具有雙曲型方程的特點。對于對流占優(yōu)問題的求解，采用常規(guī)的Galerkin有限元方法，為了避免求解結果產生數(shù)值振蕩，獲得穩(wěn)定解，則應使每個單元的局部Peclet數(shù)p此＜2

Peh=&zD這里h為單元的最大尺寸，|v|為單元中的最大速度分量值。因此，用本文方法求解對流占優(yōu)對流擴散問題，要得到穩(wěn)定解，則要通過加密有限元網格來實現(xiàn)。有限元法簡述121有限元法的發(fā)展概況血閉網有限單元法最初是在20世紀50年代作為處理結構力學問題的方法出現(xiàn)的。追溯歷史,早在1943年,庫蘭特(Courant)已應用了“單元”的概念。在1945~1955年，阿吉里斯(Argyris)等人在結構矩陣分析方面有很大發(fā)展。1956年，特勒(Turner)、克拉夫(Clough)等人把剛架位移法的解題思路推廣應用于彈性力學平面問題：他們把連續(xù)體劃分成三角形和矩形單元，單元中位移函數(shù)首先采用了近似表達式，推導單元的剛度矩陣，建立結點位移與結點力之間的單元剛度方程。i960年，克拉夫首次把這種解決彈性力學問題的方法，給予特定的名詞，稱為“有限元法在我國，60年代初期馮康等人獨立地建立了有限元方法的數(shù)學理論并應用于工程實際問題。有限單元法是基于變分原理和分片多項式插值，結合電子計算機的發(fā)展與推廣而迅速發(fā)展起來的一種求解微分方程的數(shù)值方法。由于有了現(xiàn)代化的、先進的計算工具，使得有限單元法在應用上已遠遠超過了原來的范圍。它已由彈性力學平面問題擴展到空間問題和板殼問題，能對原子能反應堆、拱壩、飛機、船體、渦輪葉片等復雜結構進行應力分析；它已由平衡問題擴展到穩(wěn)定問題和動力問題，能對結構在地震力與波浪力作用下的動力反應進行分析；它已由彈性問題擴展到彈塑性與粘彈性問題，能解決土力學、巖石力學、斷裂力學等問題；它已由結構的應力分析擴展到結構的優(yōu)化設計。除此，它在流體力學、熱傳導、磁學、建筑聲學、生物力學等等方面都有不同程度的應用。總之，兒乎所有的連續(xù)介質問題和場問題都可以用它求得滿意的數(shù)值結果?？梢灶A計，隨著現(xiàn)代力學、計算數(shù)學和計算機技術的發(fā)展，有限單元法作為一個具有鞏固理論基礎和廣泛應用效力的數(shù)值分析方法，必將在國民經濟建設和科學技術發(fā)展中發(fā)揮更大的作用，其自身亦將得到進一步的發(fā)展和完善。122有限元法解題步驟有限元法求解數(shù)學物理問題的主要步驟如下：.寫出變分表達式根據(jù)最小位能原理或虛功原理等變分原理或加權余量(方程余量與權函數(shù)正交化)原理，建立起與微分方程初邊值問題等價的變分方程(區(qū)域上的積分表達式)。這是有限元方法解題的出發(fā)點。.區(qū)域剖分為了把積分表達式離散化，根據(jù)求解區(qū)域的形狀以及實際問題的物理特點，將區(qū)域剖分成若干大小不一、兒何形狀規(guī)則的單元，并確定單元中的結點數(shù)目與位置；然后對單元、結點按一定要求進行編號，以便在單元上進行插值。.確定單元基函數(shù)根據(jù)對近似解可微性要求,選擇滿足一定插值條件的有局部支集的分片插值函數(shù)為單元基函數(shù)，確定單元中結點數(shù)目及單元形狀函數(shù)。.單元分析單元分析的目的是建立單元有限元方程。將單元中的近似解表示為單元基函數(shù)的線性組合，再將它代入積分表達式中，并對單元區(qū)域進行積分，就可獲得含有待定系數(shù)的代數(shù)方程組（或常微分方程組）。這個方程組一般稱為單元有限元方程。實際做法是采用單元形狀函數(shù)來推導單元有限元方程的系數(shù)矩陣（也稱單元剛度矩陣）和右端項（也稱為單元荷載向量）。.總體合成所謂總體合成，就是將區(qū)域中所有單元有限元方程按一定法則進行累加，形成總體有限元方程。實質是將分解開的單元積分表達式重新組合起來，形成總體區(qū)域上的積分表達式?？傮w有限元方程中的未知數(shù)正是求解函數(shù)在各個結點上的參量（函數(shù)值或其導數(shù)值）。目的就是將單元系數(shù)矩陣合成為總體系數(shù)矩陣。.邊界條件處理自然邊界條件一般在積分表達式中得到滿足。邊界條件處理主要是如何使本質邊界條件上結點的函數(shù)值滿足指定的本質邊界條件。這可以通過對總體有限元方程按一定法則進行修正而實現(xiàn)。.解總體有限元方程，計算有關物理量根據(jù)本質邊界條件修正后的總體有限元方程，是含有全部未知量的封閉方程組，可采用合適的數(shù)值計算方法求解。當求出全部待定量后，即可獲得近似解的表達式，就可根據(jù)題意計算有關物理量。面向對象程序設計簡述面向過程的編程方法和面向對象程序設計皿聞近年來，面向對象的技術正在成為一種先進的軟件開發(fā)方法，它是一種運用對象、類、繼承、封裝、聚合、消息傳遞、多態(tài)性等概念來構造系統(tǒng)的軟件開發(fā)方法，有利于克服基于過程的有限元程序在編制及維護等許多方面的不足，在系統(tǒng)的可維護、可重復使用和可擴充方面有了很大的益處，在友好界面及可視化方面也有較好的效果。目前國際上一些非常著名的有限元分析軟件大多數(shù)采用結構化程序（通常采用Fortran語言）來進行設計。傳統(tǒng)的結構化程序設計思想將整個程序分成若干個相互獨立的程序段構成一些功能模塊，每個模塊又可以由若干個子程序構成，用若干個只有一個入口和一個出口的控制結構編寫。程序的定律是：程序=（算法）+（數(shù)據(jù)結構）。即算法是程序的獨立體，數(shù)據(jù)結構也是一個獨立的整體。兩者分開設計，以算法為主。這種傳統(tǒng)的結構化編程模式被稱作是“面向過程”的編程方法。這種方法的優(yōu)點是：程序具有良好的結構，易于設計、易于理解、易于調試修改，從而提高程序設計和維護的效率。但是隨著程序編制規(guī)模的增大，其弊端便顯示出來。它把數(shù)據(jù)和過程作為兩個相互獨立的實體。在編寫程序時，必須時刻考慮所要處理的數(shù)據(jù)格式與算法函數(shù)的調用關系。對于有限元軟件開發(fā)者來說，每次增加舊版軟件的功能（如增加新的分析類型、新算法、新單元等）就要對原有程序碼進行修改和擴充。必須再次花很多精力弄清楚原有代碼和隱藏在其中的數(shù)據(jù)與函數(shù)錯綜復雜的關系，從中分辨出那些代碼在新版軟件中是可保留的，哪些是要修改和補充的，并精心對上述代碼進行重新組合，在保證數(shù)據(jù)安全調用的基礎上，盡可能提高已有代碼的重復利用率。隨著軟件功能的增強，程序代碼和編程工作量也往往大幅度增加。由于數(shù)據(jù)和代碼的分離，總存在使用錯誤的數(shù)據(jù)調用正確的程序模塊，或使用正確的數(shù)據(jù)調用錯誤的程序模塊的危險。這樣，使數(shù)據(jù)和程序保持一致便成為軟件開發(fā)者的一個沉重負擔。80年代后期出現(xiàn)的面向對象程序設計（OOP：Object-OrientedProgramming）注重于系統(tǒng)整體關系的表示和數(shù)據(jù)模型技術，即把數(shù)據(jù)結構與算法作為一個獨立功能模塊。程序定律被重新認識為：程序=（算法+數(shù)據(jù)結構）。即算法與數(shù)據(jù)結構是一個整體，算法總是離不開數(shù)據(jù)結構，算法含有對數(shù)據(jù)結構的訪問，算法只能適用于特定的數(shù)據(jù)結構。這是面向對象的程序設計的基礎，在面向對象中算法與數(shù)據(jù)結構被捆綁成一個類，從這樣的角度看問題，就不用為如何實現(xiàn)通盤的程序功能費盡心機。面向對象程序設計所擁有的數(shù)據(jù)抽象和信息隱蔽等機理以及面向對象的繼承性、封裝性、多態(tài)性等特性為軟件開發(fā)提供了理想的模塊化機制和比較理想的軟件可重用成分，現(xiàn)已廣泛應用于屏幕界面、數(shù)據(jù)庫、文件管理系統(tǒng)以及其它一些非數(shù)值計算領域。132目前國內外對OOP編程在工程計算方面的探索國外將OOP引入有限元的工作大致是從1990年開始的，它們針對有限元的典型分析步驟提出了各種各樣的類設計的方案，這些思想還不夠成熟，也不統(tǒng)一。國內在采用面向對象方法設計有限元程序方面的研究目前正在起步，僅有個別文獻涉足該領域。文獻［21］主要介紹了面向對象程序設計方法具有的意義和發(fā)展，以及面向對象在有限元程序中的應用前景。文獻［22］主要以板單元的有限元為例，闡述了采用面向對象方法和VC++語言進行有限元程序開發(fā)的過程。國內一些院校如西南交通大學、大連理工大學、同濟大學、武漢理工大學等多所高等院校也都在這方面做了比較深入和具體的工作，而求解三維對流擴散方程的面向對象有限元程序設計的文章還相對較少。本文研究的主要內容在傳熱、傳質學研究中，對流擴散問題的數(shù)值求解具有重要的理論和實際意義。目前，用來求解這一類問題的數(shù)值解法有很多，研究的角度也有所不同。對于三維定常的對流擴散方程的求解，國內外許多學者也給出了各種數(shù)值方法的求解，但是很多程序開發(fā)的具體實現(xiàn)一般都采用結構化語言，用面向對象的程序設計語言來開發(fā)相對較少。本文主要研究Peclet數(shù)較小時，三維定常的對流擴散方程的有限元解法，用面向對象的C++程序設計語言開發(fā)出有應用價值的通用程序,對相關的對流擴散問題進行數(shù)值模擬。本文的主要內容總結如下：.闡述了本文的研究背景，綜述了面向對象程序設計的基本原理。.采用四結點四面體單元和八結點六面體等參數(shù)單元進行三維定常對流擴散方程的有限元分析。.用面向對象的程序設計語言C++來具體實現(xiàn)三維定常對流擴散問題的有限元程序的開發(fā)，融入用C++編制的大型通用有限元科學和工程計算軟件。.給出了四個算例，并對分析結果進行初步討論和評價，從而驗證本文方法的有效性和實用性。2面向對象程序設計基本原理面向對象是一種試圖模仿現(xiàn)實世界的方法，它遵循認識方法學的基本原理，運用我們已形成了的對現(xiàn)實世界的概念、分類和抽象方法，將現(xiàn)實世界分成不同的對象類(class),如“人”、“汽車”、“建筑”等，不同的對象的組合及其相互作用構成了我們要研究分析和構造的客觀世界。類與對象阿2財類(class)就是種類，是具有共同行為的若干對象(object)的統(tǒng)一描述體，是對象的抽象描述。每一類對象均享有一組屬性(attribute)或行為特征，比如每個結構構件的都有長度、截面積、材料特性等。類是抽象數(shù)據(jù)類型的實現(xiàn)。一個類的所有對象都有相同的數(shù)據(jù)結構，并且共享相同的實現(xiàn)操作的代碼，而各個對象的屬性可根據(jù)施加在對象上的各種因素而變更。因此，類是所有對象的共同的行為和不同狀態(tài)的集合體。對象是一種事物，是正在創(chuàng)建的系統(tǒng)中的一個實體。從實現(xiàn)形式上講，對象是狀態(tài)和操作(或方法)的封裝體。對象對于實現(xiàn)系統(tǒng)的完整功能都是必不可少的。你永遠不會寵愛貓的定義，你只會將單個貓作為寵物。你會在貓的概念和在你的起居室到處跑的特定的貓之間化一條明確的界限。各個類的對象間通過消息進行通信。所謂消息，實際上是一個類的對象要求另一個類的對象執(zhí)行某個服務指令，指明要求哪一個對象執(zhí)行這個服務，必要時還要傳遞調用參數(shù)。系統(tǒng)功能的實現(xiàn)，就是通過一系列對象消息的傳遞，執(zhí)行--系列服務達到的。對象類和對象兩者是不同的，對象類是對具體對象的抽象(abstraction),對象是對象類的實例化(instance)。比如，你、我、他是人類，我們每一個人是人類這個對象類的單個實例。面向對象的基本特征網封裝(encapsulation)封裝是將數(shù)據(jù)和對該數(shù)據(jù)進行合法操作的函數(shù)封裝在一起作為一個類的定義。封裝的類有一個類名，程序中由類名說明對象實例。由于對象的封裝性，外界只能看到對象界面上的消息，或稱接口(功能函數(shù))，對象內部(數(shù)據(jù)和功能函數(shù)的具體代碼)對外界是隱藏的，其目的是將對象的使用者和對象的設計者分開，使用者不必知道對象功能實現(xiàn)的細節(jié)，而只需用設計者提供的消息來訪問對象。就如同你可以使用冰箱而不必關心其壓縮機的工作原理一樣，用戶也可以使用一個設計優(yōu)良的對象而不必了解其內部的數(shù)據(jù)成員。封裝的另一個好處是類的內部改造不會影響外界其他部分。繼承（inheritance）繼承提供了創(chuàng)建新類的一種方法，這種方法就是說，一個新類可以通過對已有類進行修改或擴充來滿足新類的要求。新類共享已有類的行為，而自己還具有修改的或額外添加的行為。因此，可以說繼承的本質特征是行為共享。一個類從另一個類繼承特征，稱為派生一個類，所派生的類稱為派生類或子類（subclass）。其上一?層的那個類稱為基類或父類（superclass）。類的派生過程可以無限下去。派生類從基類中繼承特征，從而實現(xiàn)已有功能模塊的可重用性。派生類在繼承基類所有方法和屬性的同時，還可以加入其自身的方法和屬性，從而實現(xiàn)功能模塊的可擴充性。比如，我們可以通過建立一個單元類（基類）來派生出各種具體的單元類（子類），如平面單元類、梁單元類等等，而平面單元類（這時是基類）又可以派生出三結點平面單元類、四結點平面單元類等子類。繼承可以讓你重用父類的代碼，專注于為子類編寫新代碼。多態(tài)性（polymorphism）多態(tài)性就是同種操作具有多種形態(tài)，在C++中，多態(tài)性乂被直觀的稱為“?個名字，多個函數(shù)它主要是通過函數(shù)名重載實現(xiàn)的。重載函數(shù)是指程序中對同名函數(shù)進行調用時，編譯器會根據(jù)函數(shù)參數(shù)的類型和個數(shù)，決定該調用哪一段函數(shù)代碼來處理這個函數(shù)調用。在C++中，多態(tài)性的實現(xiàn)與聯(lián)編這一概念相關。一個操作函數(shù)的調用與適當函數(shù)代碼相對應的動作，叫做聯(lián)編。聯(lián)編分為靜態(tài)聯(lián)編和動態(tài)聯(lián)編。在程序編譯時就確定哪個重載函數(shù)被調用的，稱為靜態(tài)聯(lián)編。直到程序運行時才能確定哪個重載函數(shù)被調用的，稱為動態(tài)聯(lián)編。多態(tài)性也分為靜態(tài)和動態(tài)兩種。靜態(tài)多態(tài)性是指定義在一個類或一個函數(shù)中存在的同名函數(shù)，它們可根據(jù)參數(shù)表（類型及個數(shù)）區(qū)別語義，并通過靜態(tài)聯(lián)編實現(xiàn)。例如，常用的加法，整數(shù)相加得到整數(shù)，復數(shù)相加得到復數(shù)，矩陣相加得到矩陣。同樣是加法操作，針對不同的對象可有不同的操作形態(tài)而產生不同的執(zhí)行結果。動態(tài)多態(tài)性是指定義在同一個類層次中的不同類中的重載函數(shù)，它們一般具有相同的參數(shù)表，因而要根據(jù)指針指向的對象所在類來區(qū)別語義，并通過動態(tài)聯(lián)編實現(xiàn)。例如，對于從同一個基類圖形類派生出來的兩個子類三角形類和圓形類，有名稱相同的操作（函數(shù)）求面積，但三角形和圓形要采用各自的計算公式。10面向對象的程序設計?算法與數(shù)據(jù)結構面向對象的程序設計的本質時把數(shù)據(jù)和處理數(shù)據(jù)的過程當成一個整體——對象，因此，面向對象的程序設計的基礎是類。類是C++的基本單元，它是由結構型數(shù)據(jù)演化而來的。結構型數(shù)據(jù)是將事物用兒種基本數(shù)據(jù)類型作為整體來描述某一事物。如人的結構型數(shù)據(jù)可描述為：structperson{charname[20];intage;floatweight;floathight;卜結構型數(shù)據(jù)是一般結構化語言對事物的描述方式，它比用基本的實數(shù)、整數(shù)來描述事物要先進一步。類是結構型數(shù)據(jù)的進一步演化。類結構是將描述事物的數(shù)據(jù)和處理數(shù)據(jù)的函數(shù)作為整體來描述事物。圖2.1算法與數(shù)據(jù)結構關系類結構是軟件方法的發(fā)展必然。較早的軟件開發(fā)用結構圖2.1算法與數(shù)據(jù)結構關系程序=（算法）+（數(shù)據(jù)結構）。即算法是一個獨立的整體，數(shù)據(jù)結構也是一個獨立的整體。兩者分開設計，以算法（函數(shù)或過程）為主。見圖2.1o圖中的虛線表示上面的算法與下面的數(shù)據(jù)分離。雙箭頭表示作為黑盒的函數(shù)輸入輸出的數(shù)據(jù)。隨著時間的流逝，軟件工程師越來越注重于系統(tǒng)整體關系的表示和數(shù)據(jù)模型技術，即把數(shù)據(jù)結構與算法作為一個獨立功能模塊。程序定律被重新認識為：程序=（算法+數(shù)據(jù)結構）。即算法與數(shù)據(jù)結構是一個整體，算法總是離不開數(shù)據(jù)結構，算法含有對數(shù)據(jù)結構的訪問，算法只能適用于特定的數(shù)據(jù)結構。這是面向對象的程序設計的基礎，在面向對象中算法與數(shù)據(jù)結構被捆綁成一個類，從這樣的角度看問題，就不用為如何實現(xiàn)通盤的程圖2.2構成程序的對象Fig.2.2Theobjectsofprogram序功能費盡心機。現(xiàn)實世界本身就是一個對象的世界，任何對象都具有一定的屬性與操作，也就總是能用數(shù)據(jù)結構與算法兩者合一地來描述。這時候，程序被再次另眼相看：對象=(算法+數(shù)據(jù)結構)，程序=(對象+對象+…)。即程序就是許多的對象的組合，而對象又是一個個程序實體。人們不再靜止的去看待數(shù)據(jù)，而把它看成是程序的單位，一個程序的分子或對象。它本身不但包含算法而且包含數(shù)據(jù)結構。對象中可以包含小的對象，它也可能含于其它對象之中。見圖2.2o2.3.2面向對象分析和設計.面向對象分析(OOA-objectorientedAnalysis)面向對象分析的目的是完成對問題空間的分析和建立系統(tǒng)模型。其具體任務是確定和描述系統(tǒng)中的對象、對象的靜態(tài)特性和動態(tài)特性。對象間的關系及對象的行為約束等。其主要內容歸結為靜態(tài)結構分析和動態(tài)行為分析。a.靜態(tài)結構分析靜態(tài)結構描述對象、類以及類間的相互關系。靜態(tài)結構分析首先要確定對象和對象類。對象是一種具有簡明界面及應用意義的實體的抽象。所有對象都具有唯一的標識，可予以區(qū)分。對象類描述的是一組具有共同屬性、操作和語義特征的對象，其中屬性是對象的性質，操作是對象的行為。在確定了對象和對象類的基礎上，還要確定對象類之間的關系。對象類之間的關系主要有：層次關系和聚合關系。在分析過程中根據(jù)類的共性及個性將類組織成不同層次。高層次的類表達共性形成基類；低層次的類表達個性，形成子類。子類通過繼承機制來獲得基類的屬性和操作。這就是層次關系。聚合關系是對象之間的組合構造關系。例如，單元對象中涉及到剛度矩陣的計算，同時單元對象的許多數(shù)據(jù)都來自結點對象，因此單元類中就包含有矩陣類和結點類，當然還有材料類、荷載類等。b.動態(tài)行為分析動態(tài)行為描述了系統(tǒng)中對象的合法狀態(tài)序列，它包括兩方面內容：一是單個對象自身的生命周期演化；二是整個對象系統(tǒng)中對象間的消息傳遞和協(xié)同工作。這部分內容需要結合具體程序說明，此處不再詳述。.面向對象設計(OOD-ObjectOrientedDesign)面向對象的設計是確定問題的解決方案的過程。面向對象設計的涉及面有大有小，結構分析軟件的設計涉及用戶界面、系統(tǒng)數(shù)據(jù)結構、系統(tǒng)資源分配、數(shù)據(jù)交換存儲、圖形顯示等許多方面。綜合起來，面向對象設計包括對象設計和系統(tǒng)設計兩方面的內容。a.對象設計12對象設計以對象分析為基礎，分析過程已給出對象模型。在設計對象時可根據(jù)需要進行擴展，增加和優(yōu)化屬性數(shù)據(jù)及功能操作，對分析中為考慮周全的類關系進行調整，盡可能利用繼承的優(yōu)點，提高代碼的復用程度。b.系統(tǒng)設計系統(tǒng)設計是為了實現(xiàn)需求目標而對程序的系統(tǒng)結構進行的總體設計。通過程序系統(tǒng)的調用各對象功能操作，達到解答問題的目標。通常系統(tǒng)結構本身也包含在對象內部。系統(tǒng)設計主要包括下列內容：數(shù)據(jù)的出去和處理設計，用戶接口及相應過程設計，系統(tǒng)組織設計等。2.4面向對象程序編程和開發(fā)工具面向對象的程序設計語言有許多，其中的C++以其高效實用性被應用于程序設計的眾多領域，功能強大的開發(fā)環(huán)境和許多技術成熟的應用類庫使C++編程效率大大提高。許多程序的開發(fā)和已有程序的二次開發(fā)都采用了C++。比如，Windows程序的編寫，AutoCAD二次開發(fā)工具ObjectARX等。C++由C改進擴充而來，包含了C的全部特征，添加了對面向對象編程的全面支持。主要包括：可進行類定義；為類提供了數(shù)據(jù)操作的函數(shù)及操作符；為類提供了構造方法及析構方法；對類內部的數(shù)據(jù)結構進行了封裝性的劃分；提供了對繼承關系及多重繼承的支持；提供了對多態(tài)性及函數(shù)重載的支持等。因此，本文采用VisualC++6.0來進行三維對流擴散問題的有限元分析程序的開發(fā)。133三維對流擴散問題的有限元分析對流擴散現(xiàn)象在研究環(huán)境保護的流體力學中是經常遇到的。流體中的污染物包括化學的、生物的或熱的污染，在流體介質中的彌散是環(huán)境保護問題中越來越重要的研究課題。本章將研究三維定常的對流擴散方程的有限元數(shù)值解法，主要闡述三維定常的對流擴散方程的有限元分析過程。概述三維對流擴散方程對流擴散方程是流體力學中的典型方程，連續(xù)方程、動量方程、能量方程都可以寫成對流擴散方程的形式。根據(jù)方程(1.1),在定常的粘性不可壓縮流體流動的對流擴散問題中，假定求溶解于流體中的物質的濃度U在流動過程中的變化規(guī)律,則U所滿足的控制方程為：TOC\o"1-5"\h\zJ(包、 8uD ——ku=finQ,i=1,2,3, (3.1)&Idxi) 加邊界條件為：U=UoonTi,(3.2)duq=Di ni=q。ox,on「2,(3.3)TOC\o"1-5"\h\z其中具有光滑邊界,r.ir2=0的區(qū)域CU33R3,邊界「上的單位外法向r r，幾，n ,v,v數(shù),物質的分解速率系數(shù)為k,這里〃。，q。，/為已知函數(shù)。本文不妨設方程(3.1)中的坐標系(Xi,X2,X3)為三維笛卡兒直角坐標系(”z),則在直角坐標系(蒼乂z)中，三維定常的對流擴散方程為：d(du、d( 3( duonTi, (3.5)onronTi, (3.5)onr2, (3.6)u(x,y,z)=uo(x,y,z)du a“ a”p.dx dy 6z聯(lián)立(3.4)、(3.5)和(3.6)式，就是三維對流擴散問題的定解問題。此時具有光滑邊I-Vi量〃=(〃l23),流場中流體的速度為V=(V.23),物質在三個方向上的擴散系I-ku=/inQ,I-ku=/inQ,(3.4)山+。 盟+。*妊伏(工,y,z),v,v ,D,D質的分解速率系數(shù)為&,且有〃。，歟，/為已知函數(shù)。3.1.2變分方程.散度定理設「是m維空間區(qū)域C的邊界，且邊界「是分片光滑閉曲面。如果函數(shù)P(X)[i=l,2,L,zn,X=(孫X2,Lx“)]及其一階偏導數(shù)在閉區(qū)域。上連續(xù)，則有m mcos(n, (3.7)則可寫成向量形式(3.8)此式稱為散度定理(Divergencetheorem).變分方程設任意口(x,y,z)wCi(C),且(p(x,y,z)h=00用中乘方程(3.4)的兩端，并在C上求積分。)(pdfldm d)(pdfl(D，dxi dyi dZ2這里，d2Ud2iid2U0x2dyi0Z2dd(u此時根據(jù)散度定理(3.8),則可得到，I(Dm,n.+Du,小+°必)峭r界「=「11￡,rjr2=0的區(qū)域qu*R3,邊界「上的單位外法向量〃=（%,：）,流場中流體的速度為「=（口>「），物質在三個方向上的擴散系數(shù)（D,,力，物0普dC=j￡p，）d「n是「的外法線。式（3.7）是m維情況下的高斯（Gauss）公式。如果用P表示向量,VgPJQ=\pgndV?▽”?(pdC-JktupJQ=J/(pdQ+D\u?▽〃?(pdC—Jku中dQ=J/(pdC+Dy+Dz豆3<p）+用軀。-J(-J(Dxxxyyy而q)(x,y,z)h=0,所以1q(pd==1q(pd「+Lq(pd「=Ji：qo(pd「因此，邊值問題(3.413.6)的變分方程為，(￡)vWt(Pv+￡)v?y(Py4-D〃傘)dfl(3.9)四面體單元為了把變分表達式離散化，根據(jù)求解區(qū)域的形狀以及實際問題的物理特點，將區(qū)域剖分成若干大小不--、兒何形狀規(guī)則的單元，并確定單元中的結點數(shù)目與位置；然后對單元、結點按一定要求進行編號，以便在單元上進行插值。首先對三維對流擴散問題采用四結點四面體單元進行單元分析，建立單元剛度矩陣與單元載荷向量。321位移函數(shù)最簡單的三維單元是四結點四面體單元。我們對區(qū)域c作四面體單元剖分，單元e.(”=1,2,L,NE),結點(四面體頂點)p,(x,,y,z，)(i=1,2,L,NP)。設m(X/,y-,Zi)=w(pi)=m,(i=l,2,L,NP),任取單元e，它的四個頂點是pt,P),pm,pi,記為e=(pi,pj,pm,pl)o它們的排列順序按右手準則，即p>,pj,pm按右手螺旋旋轉時，p,在拇指方向(如圖3.1所示)。P'P,J(D邛甲武?！暧?方8也?|\，?Vw,(pdQ-Jku^dQ=J/(pdC^Dxxxyyyzzz^Dxxxyyyzzz)</Q+pmv-V”?(pdC+Jku(pdQ=Jq(pdf-]/(pdC圖3.1圖3.1四面體單元Figure3.1Thetetrahedralelement由u在四頂點的值，恰好祚單元內可以隼定一線性單值。即茸令

+Jy?V〃(pdQ+Jk〃(pdQ=Jq()(pdT-jf(pdQu=aix+a?y+oi3z+oi4 (J.1U)16則它滿足[m；=CLtXi+CLiy,+a.3Zi+cu\uj=aixj+a.2yj+a?zj+oi41ui=aix/+azy/+a?z/+ou解出ai、a?、ou、a4代入(3.10)再適當整理可得U=Ni〃i+NjUj+NmUm+NlUl=IN]{3}e (3.11)其中1xyz 1xyt11Xj yjh 11x yINi= ， Nj= ，6Ve1Xin y(n& 6Ve1Xm Z?ix， y ： i x， y ：1Xi yt v 1 Xi y> Zi11Xj yj: 11Xj yjZjNm= ， Nl=6Ve1X y' 6Ve1Xm Z>nt\Xi yi 1 yzX此處匕代表四面體e的體積，1Xi y. Zi11Xj Vj ZjV,二一 (3.12)61Xm ym Z"i1Xi yi Zt[N]=[N Nj iV.Ni]{8}e=[Ui Uj Um Ul]r若記| 「 iM=6],[ax+b，y+c，z+d<] = (3.13)則1 yj Zj 1 Xj zjUi——1y>n Z>n9 bi—1 XmZ91 yi Zl 1 Xi zi1 Xj yj Xj yj ZjCi——1Xm ym, di—Xm y”i1 Xi yi Xi y, N其余的G,b,,a,乩(S=_/,"/)亦可照此排列方式類似寫出。但必須注意不能I=aiXm+a2 +ouzm+ou

Clmdlj像平面問題那樣循環(huán)轉，因為它不能保證適合右手準則。例如就不適合。凡在原排列中排序為偶數(shù)的(例如j排序為2),若以它為首開始循環(huán)，則剛好反向，這時只要按原公式乘(-1)Clmdlj1ymZm1XmZmaj=1 yiZi，bj=-IXiZi，1yt尊1XiZi1Xm*XmymZfMCj=1 Xiy/，dj=Xi1日,1XiXy霹Zi1XlZIam=1y>Zi，bm—1XiZi，Zi1XjZj1XiVXIyiztCm=—1Xiy，dm-Xi■Zi,1XjyiXjZj1yZi1XiZiai=\ yjZj，bi=-1XjZj，1ymZm1XmZm1XiyXi小Cl=1 Xjy，，di=XjyjZJ1XmymXmyM的梯度向量V"=[g"dxduSy]可表為：「1fjMSNidNmM1Idxdx11SyduI1dNidNjdNmSydy洲1{5b=[B]{5}e1(3.14)1. dN_idNmI1dzdzi不妨令「dN,LdeGNidedNmdNi^de(c=X,y,z)(3.15)則叼=6；依durdzTOC\o"1-5"\h\zII\dx dxI尸產 ay I_j迪 dNi IIL&U\_dz dz JdeJ

TOC\o"1-5"\h\z\o"CurrentDocument"1 I[^]=e[b.bjb,nbi], [B：]=[aCjCmct]ol> OK「BJ b|「8、[「加[fi]=11By111 UILb.-U因此，■:

I1118-31118-3一〃一(3.16)一=[B.]{3}<.(c=x,y,z)曲(3.16)322單元剛度矩陣變分方程(3.9)按剖分單元改寫為NE NEu(p+￡)〃(p+?！?pen en(3.17)單兀e：ptp>p?,pi,函數(shù)“(x,y,z),別為Us> ?,?,o〃(x,y,z)=[N]{3}‘(p(x,y,z)=[AT]{8?}?r”r[B-]{8}ju1鬧1111l}11Ldz1JJ(Dm(P. 四 (9>?+D(9)tii+D(仰itT T Tef* 3 5rr rCU CU1 ；+V，7(p(x,y,z)在結點p5(s=i,j,m,l)上的函數(shù)值分{8}e=[UiUjUmUi]t*Um*u*i]rrd(p11「田,]{3小1,L1*1SyIL5zIJ*fl+v；,)dxdydz2J[Dxxxyyyn=l其中rnIdxI|&|=|B|{5}=|[fi]{5}IILr-UIL[r-"」zz二)dC+ZJv-Vw-(pr/Q，=1+ZJku(pdQ=ZJ<7?<pdr-ZJ/(pJQn=l n=l“I「2,若「2與 不相交，則jqo(pd「=0{3*}f=[uiu*IdxI IdxIlal=l[8]{3}I V(p=| |=|[Bv]{5},IIL&'U 曳|I[叼{&J」+Oyyy工工二)dxdydz~j(。*xxyyyzzz)(lxdydz=@}〃[J(dM]W丹2【4][p[+2[8][8])dxdydZ]i22s2oyv-Vw?(pdxdydz=TOC\o"1-5"\h\z={3，}小],{3}‘這里[k]e為單元剛度矩陣T T Te(3.18)T T T T(3.18)e e很顯然伙卜是4x4矩陣kitekije kime hIK kejj kejm,kmjeLkli履h||4x4h=I(D，dxdx0ydy

dzdz)dxdydz+vN323單元荷載向量f<pdxdydz=I(<p)rfdxdydz={6>}tlf[N]rdxdydz={5-}r{F},(3.20)(3.21){Fi}r=I/[]rdxdydz=[Fi,F\e)F\(mFnr(3.20)(3.21)Fi?=[fN,dxdydz=如果單元e不是邊界單元或者單元e是邊界單元，且若「與de不相交，則單元荷載向量{尸}，=-{尸?。喝绻麊卧猠是邊界單元，且若R與灰相交，設r=de\Vi是其邊界，則有o qods={3,}rJ飲[N]rds={3*}r{B}c={3.}〃[J(v』N]叫「[NJ[BJ+“N][B])dxdydz\{5}Jk〃(pdxdydz=jk(<p)udxdydz(Ll[N][N]dxdydz){6}，[(DvK<P即P〃*10(M<PVw?(p)+(Z〃(p)]dxdydz={6?卜[][B.]+Dt[B>][]+D[B；][B；])dxdydz+](v.[N][B,]+v,[/V][B,]+v；[N][B：])dxdydz+\k[N][N]dxdydz]{^>}，=J(dI4][41+0UJ[B]+D[微分以+J(v,[N]x[B]+v[/V][B]+J]^^|此+\k[N][N]dxdydzkmlc|單元剛度系數(shù)k為如辿+於@」加迎f8N,dN,dN, f+JW- +丫… ：， )dxdydz+JkN,,dxdydz(s,t=i,j,mJ)(3.19)[q(pd「=J(<p)[^]rJ5=[FF勿 [^]rJ5=[FF勿 Fie1]r(3.22)"Mze(3.23)Fze=JrOoN,ds(s=i,j,m,l)如果單元荷載向量統(tǒng)一記為｛用‘，則有(3.23)(1)如果單元e不是邊界單元或單元e是邊界單元，且若「與良不相交，則｛尸一｛E｝,；(2)如果單元e是邊界單元，且若口與加相交，貝IJ｛用,=｛氏｝‘-｛R｝,。整體分析為了形成總剛度矩陣和總荷載向量，把單元剛度矩陣，單元荷載向量的表達式代入(3.17)式，為了便于疊加，對｛3｝,,｛用，和［kh進行擴充，擴充為NP維向量和NPxNP階方陣。對某一個單元e來說，頂點p,,p/,p“,p,保持為原如圖順序，但編號也可能不依次了，即i,/大小不定。如j<i<m<I其中沒寫完整，中間其它不是系數(shù)的為Oo在擴充以后，我們可以疊加：(3.24)NE NE〃=l n=l元剛度矩陣和單元載荷向量的元素按腳標編號對號入座放到總剛度矩陣和總載荷向量的相應位置上去，進行疊加而成。leiIl 的l11TOC\o"1-5"\h\z{用,=|M| [口‘="’ ’e eU LLJ￡{>}〃"[&]“{3},"=Z{3}“”{F},“加}(Z伙卜){3}={3?}(Z{用 Q{5-}([*]{5}-{F})=0這里伙］=2伙］“ {F}=Z{F}，“分別為總剛度矩陣和總荷載向量。它們都是由單約束條件的處理自然邊界條件一般已在變分方程中得到滿足，本質邊界條件需要通過對總體有限元方程的修正，使其滿足。實質上是在近似解表示為總體基函數(shù)線性組合的表達式中構造一個滿足本質邊界條件的特解。如果結點p,G「(r=1,2,L,NP)即要使u(pr)=Uo(pr)@UOr則 l：=k,r0:=的，& (_/=l,2,L,NP&#r)r

QF-kirUo:=F (i*r)如果另有點，繼續(xù)采用類似的方法來處理。按此想法，約束邊界條件也可以這樣來處理，就是將已給出本質邊界值UOr的所在行列的全部元素去掉，右端項元素中將第r個元素去掉，其余元素仍按前面方法處理為F-k,rU0(i=l,2,L,NP,且,如果有兩個以上本質邊界結點，則按上述規(guī)則逐個進行。這個方法可以使系數(shù)矩陣的階數(shù)減小，但是由于縮減后的矩陣，原先的下標序號已經打亂，應重新編號。最后，解線性代數(shù)方程組，求結點的值。再求單元內u(x,y,z)的值?(x,y,z)=Nm+Njuj+L+Niui (x,y,z)we“，但是u|n=uo。八結點六面體等參數(shù)單元對于空間問題，經常使用等參數(shù)單元，而且效果顯著。這里研究經常采用的八結點六面體等參數(shù)單元。從局部坐標系化E")下的母單元開始，它是一個將單元的中心設為坐標原點、邊長為2的正立方體，正立方體的八個頂點取為單元的結點，如圖3.2所示。Iq 1m n圖3.2六面體單元Figure3.2ThehexahedronelementU：-Fru(x,y,z)=[N]4mawk2wf (3.25)其中，[N]=[N,NiNtMNmN?NPTV,]{8}<?=[ttiUjUkUiUrnUnUpUi]r并且，N*=N.也,T|,q)(s=i,j,k,l,m,n,p,q),M(3,T], (s,t=i,j,k,1,m,n,p,q)10(5豐t)1(s=i,j,n,p,q)(3.26)8在空間任給八個對應結點的整體坐標值,則由上述形函數(shù)便可決定一個八結點等參數(shù)單元。相應的坐標變換是1 (3.27)I其中(x,,ys,z.)(s=i,j,k,l,m,n,p,q)是給定的單元八個結點的整體坐標。當一個自變量固定時，四位移函數(shù)是另外兩個自變量的雙線性函數(shù)，因此，在立方體單元的每個側面上，該插值函數(shù)由該側面上四個結點的函數(shù)值惟--確定。因而，在相鄰單元的公共面上，只要在其四個公共點上的函數(shù)值相同，位移函數(shù)一定是連續(xù)的，即在局部坐標下滿足協(xié)調性的要求。由插值函數(shù)式的協(xié)調性即可保證坐標變換式的協(xié)調性，同時保證了整體坐標下位移函數(shù)的協(xié)調性。?的梯度向量V“=[乎]可表為：oxdy「色1「磔dNidNkdNidiLndNndNj,碼]18xdxdxdxdxdx1du11dNi8NidNkdNidNn,dNndNP碼'/XISyBySySySySy118N,dNkdNidN.dNn迦1.1dzdzdzdzdzdz1=[B]{8h一“ dN,迦dNkdN1SNadNndNr\遠不妨令=[_ _ _ —oc de oc dededededeJ(c=x,y,z)(3.28)m&,rj,q)=(i+)(i+H'H)(1+GC)fx=Zm化,n,q)x,[)，=ZM,《,n，4)y:尿ji；i：dxI

dydxI

dyI&JIdxI|qxIcbII?'IL^UL-

r8j網=11By||ILb-Ub1「8Ju因此，I1du11Ldz\一du~=[Be]{8}e(c=x,y,z) (3.29)在空間問題中，位移u是整體坐標x,y,z的函數(shù)，又x,y,z又是局部坐標匕,r),。的函數(shù)，形函數(shù)N,化,“工)是用局部坐標給出的，根據(jù)偏微分法則，可知(3.30)同理可寫出皿制如，集合起來，得到3r| /「SNsl「@「SNsl「@I11IdMII_6xIII優(yōu)出dx\dzl11 1II加11 1111dN,I1 1UII出dz11IIdNs1111加1%II&11dz1(3.31)由于形函數(shù)是用局部坐標給出的，上式左端可根據(jù)Mg,n,c）求出。另外，由坐標變換公式可求出矩陣叫，此矩陣稱為雅可比矩陣，即[8x出&~l「一皿1小1她 型\dx西dy奧11dz11-I *加la心恥彷｝=囑］㈤118r\dv\西111b：11*M｛弭IJl@出dz\\_dNs dN, 洲,1叉照」L-424dN,dN，dxdNsdydM&閨優(yōu)I西1=1的La;J[%&~1「8Ml 「閨優(yōu)I西1=1的La;J[%dE,Idx?-dx-\Sr|IIII""分?-JLJLJ1死11Z比X5X比X比[小11=1Z的■工即Xsr|L-IlZ/Xs工式Eaqj「一狙T膽11 比LIyi必1I 狙L制為1L !1yjziLL|1=M[x](3.32),狙1L丫口411生或生IJrdN,dNi皿1|更奧IIiytz11dN迦1dxy為L 1，1.t zLLII(3.33),迦T IhzILa;買Lku對矩陣［J］求逆后，由式(3.31)可得到形函數(shù)在整體坐標中的導數(shù)如下：「Ml1 1「迦]1 11 11dN.|=1 11 1J-JdNs|1 1(3.34)111dN?11 11年I根據(jù)(3.17)式，與四面體單元類似，可求出采用八結點六面體等參數(shù)單元時的單元剛度矩陣和單元荷載向量。1比1則1Sr]比11「的dN,IIXj3r)||LdNdN,IILx? J1%[y]=i(Ml y」||劭I葉心Ml 必,UIdyI[]Mt]IILazUILazUL區(qū)」4面向對象的三維對流擴散問題有限元分析程序實現(xiàn)面向對象程序設計方法是一種新興的程序設計方法，或者說是一種新的程序設計范型。該方法強調直接以問題域中的事物為中心來思考問題、認識問題，并根據(jù)這些事物的本質特征，把它們表示為系統(tǒng)中對象，作為系統(tǒng)的基本構成單位。它對描述事物的屬性、特點及事物之間的關系具有一套較完整的系統(tǒng)分析方法，并填平了語言上的鴻溝。由于本論文是在袁政強老師所編的有限元計算程序平臺上進行的，其程序類已在口7][23]中有較詳細的介紹。本章簡要介紹數(shù)值計算類、數(shù)據(jù)結構類和有限元分析程序的基本框架，系統(tǒng)地闡述采用面向對象的程序設計語言C++實現(xiàn)三維對流擴散問題有限元程序的開發(fā)過程。數(shù)值計算類由于在有限單元法計算中要大量處理應力張量、應變張量、向量、矩陣等數(shù)據(jù)。根據(jù)需要，在編制有限元程序時，編制了通用的數(shù)值計算類，包括向量類(Vector),矩陣類(Matrix)、點向量類(Point)以及張量類(Tensor)等。這使得對張量、向量、矩陣、的運算就像我們通常習慣的變量那樣簡單明了，同時使程序結構緊湊、代碼簡潔、可讀性強。同時，這些類的定義大都采用模板類，這樣可以勝任各種類型的數(shù)據(jù)，如int,float,double等，程序可根據(jù)要求自動生成相應數(shù)據(jù)類型的向量、矩陣及張量等對象，從而使程序的通用性更好。數(shù)據(jù)結構類利用數(shù)組或順序方式來組織數(shù)據(jù)結構，其優(yōu)點是存儲利用率高，存儲速度快。但是，對于數(shù)據(jù)元素個數(shù)動態(tài)增長的情形，由于數(shù)組的元素個數(shù)不能自由擴充(除動態(tài)數(shù)組外)，一旦空間用完就不能再向里加入新的元素，否則就會導致系統(tǒng)停工。另外，在表中插入新元素或刪除無用的元素時，為了保持其它元素的相對次序不變，平均需要移動一半元素，效率很低。對于有n個長度變化的有序表分配到同一存儲空間情況，按每個表的最大可能長度分配空間就會造成空間的巨大浪費。為了解決這些問題，程序采用了鏈表(List)類，鏈表結點類(ListNode)以及鏈表游標類(Listiterator)?；境绦蚩蚣芨鶕?jù)結構有限元分析方法，在具體結構分析過程中，可以將整個結構作為一個對象，稱作整體結構對象，它包含了多種成分對象。將成分對象抽象歸類時，26要遵循各類“分工明確，各行其責''的原則，明確對象類所描述的內容及其操作功能。根據(jù)這一原則，從整體結構中可細分出單元對象、結點對象、荷載對象和材料對象等。對這些對象進行描述和建模，形成所謂的結構對象類。結構對象類由抽象單元類(Element),結點類(NodeData),荷載類(Load),抽象材料類(Material)和有限元方法類(FinMethod)構成。(1)有限元方法類有限元方法類是有限元分析最基礎的類，它提供了有限元分析所需的最基本的內容：分析數(shù)據(jù)和分析方法。因此，有限元方法類的成員變量就是結點數(shù)據(jù)、單元數(shù)據(jù)、材料數(shù)據(jù)、荷載數(shù)據(jù)、約束條件等，它們分別由下面介紹的各個類提供。大多數(shù)以鏈表的方式存儲在有限元方法類中。有限元方法類成員函數(shù)則是有限元分析最基本的一系列步驟，如數(shù)據(jù)輸入、計算單元參數(shù)、計算單元剛度矩陣、單元剛度矩陣進入總剛度矩陣、計算等效結點力、計算方程組右端力向量、解方程得到結點位移、結果后處理等。結點指針(1)指針(2)指針(3)指針(4)指針(5)指針(6)指針(n)指針NULL單元(1)(2)(3)(4)(5)(6)(n)NULL指針指針指針指針指針指針指針指針材料(1)(2)(3)(4)(5)(6)(n)NULL指針指針指針指針指針指針指針指針圖4.1結點集合、單元集合、材料集合之間的關系圖Fig4.1Graphofarelationwithnodesset、elementssetandmaterialset有限元方法類中的結點集合類是有限元單元的兒何結點的集合，包括了分析模型中的所有結點，以鏈表形式存儲，稱為結點總鏈表。而每一個單元中都包含了一個存儲本單元所有結點的鏈表，稱為單元結點鏈表，它通過初始數(shù)據(jù)(結點編號)與結點總鏈表中相應結點建立一一對應的關系。因為鏈表中存儲的都是結點指針，所以這種對應關系實際上是共用同一段內存地址。這樣，無論在后面的有限元分析中結點總鏈表發(fā)生插入或刪除操作，因而改變結點的編號，每個單元都仍然可以根據(jù)結點指針準確的在結點總鏈表中找到自己的結點。這種對應關系見圖4.1o因此，這是一種動態(tài)對應關系。同樣，每一單元還對應一個單元所使用的材料類，也是通過與結點相似的方式建立的。方法類中的函數(shù)只能對該類中的數(shù)據(jù)進行操作，不能對它數(shù)據(jù)的下一級數(shù)據(jù)進行操作，這種層次性的數(shù)據(jù)操作體現(xiàn)了C++對數(shù)據(jù)的先進管理方式。如方法類中的總體剛度矩陣計算，只對單元集合鏈表進行操作。如果統(tǒng)一對單元類發(fā)出計算單元剛度矩陣并放入總剛度矩陣中的指令：’lei.Next()->GetStiffMatrix(K,*F),,單元類就會根據(jù)單元的類型，采用相應的單元剛度矩陣計算函數(shù)進行計算，計算出的單元剛度矩陣隨著單元類型的不同有不同的階數(shù)，然后根據(jù)單元所對應的結點放入到總剛度矩陣之中。在方法類中并不關心單元的類型是什么以及采用何種方法執(zhí)行指令，這些都屬于它的下屬單元類數(shù)據(jù)處理的事。方法類中采用這樣的處理方式使程序的編寫和調試都變得簡單容易，且函數(shù)的結構形式也完全相同。(2)抽象單元類抽象單元類包含了具體單元所共有的數(shù)據(jù)、函數(shù)和函數(shù)接口。其任務就是派生具體的單元類以實現(xiàn)C++的多態(tài)性，以及在有限元方法類和具體單元類之間建立聯(lián)系橋梁。首先，具體單元有不同的類型，具體類型的單元都是抽象單元的繼承。C++的繼承性允許把對象作為它本身的類型或它的基類類型來處理。這就意味著允許很多類型(從同一個基類派生的)被等價地看待，就像是一個類型，而且允許同一段代碼工作在所有這些不同類型上。虛函數(shù)(函數(shù)接口)就是這樣的一段程序，而這個基類派生的不同具體類型之間的區(qū)別，就在于每個具體類型中都含有這個虛函數(shù)的一個具體重載函數(shù)，只是函數(shù)的實現(xiàn)不同，這樣，C++的多態(tài)性會根據(jù)具體的對象類型調用相應的重載函數(shù)，這種編程方式更符合我們的邏輯思維習慣。如虛函數(shù)Stiff()剛度矩陣的計算，對于不同的單元，剛度矩陣的計算方法和計算的結果都是不相同的。但在單元集合類中并不區(qū)分單元的類型，都按統(tǒng)一的方式混合放入單元集合中。當調用剛度矩陣計算時，C++的多態(tài)性會根據(jù)具體的單元調用具體的單元剛度矩陣計算函數(shù)。虛函數(shù)在抽象單元中是不用編寫的，它只是具體單元的計算接口。在編寫具體單元時，才對虛函數(shù)進行重載，不同的重載就是不同單元的計算方法。在抽象單元中給出了許多這樣的虛函數(shù)計算接口，如：單元剛度矩陣計算等。其次，通過鏈表形式來處理單元，可以方便的實現(xiàn)單元在計算過程中的增加或使單元退出計算。在計算過程中，有些單元會因為壓碎失去剛度而退出計算，有時會因為單元裂變由一個大單元裂變?yōu)閹讉€小單元，單元會在計算過程中不斷的增加。由于鏈表結構方便的增加減少方式為處理這類問題提供了極大的方便。下面介紹抽象單元類的成員函數(shù)和成員變量。抽象單元類的成員函數(shù)除了一些輔助功能函數(shù)外大多是虛函數(shù)，它們是具體單元的計算接口，都要在具體單元中重載。由此可以看出，抽象單元類是具體單28元類和方法類之間的橋梁，通過它來實現(xiàn)C++的多態(tài)性和繼承功能。采用從基類按層次繼承來建立單元類的方法，除了可重用父類的代碼以外，另一個很重要的優(yōu)點是可以運用面向對象計算中的多態(tài)性，即動態(tài)聯(lián)編技術。通過指向單元基類的對象的指針，調用各派生單元類的成員函數(shù)。這樣做的好處在于，需要增加新的單元類時，只要在其他相關類中增加新單元類的說明，即包含類定義頭文件，而不必改動相關的類，從而大大提高編程效率及代碼的重用率。(3)結點類結點類封裝了結點的信息，包括自由度、坐標、結點位移、約束條件等。結點序號是在輸入結點時的排列順序，由于在計算過程中要不斷加入單元和結點，因此不能直接依據(jù)這個序號來將單元剛度矩陣疊加到總剛度矩陣中去。為了解決這個問題，引入了結點在方程中的開始編號Nequation,由程序自動計算，它的作用就是指示結點的每一個沒有受約束的自由度正確的進入總剛度矩陣。每一個結點存儲著自己的約束條件，分為邊界條件標識和邊界條件值。有約束的自由度標識為1,否則為其它，在有約束的自由度上給出約束值。(4)抽象材料類每一種具體單元在應用中都需要賦予一個具體的材料屬性，這用材料類來表示，正如上面介紹的抽象單元類中含有一個材料類指針，就是用來實現(xiàn)這個屬性。材料的種類繁多，其中有些材料又有不同的等級，可以細分成不同的材料。既便是同種材料，在分析過程中也可能發(fā)生變化，如材料非線性問題中的材料退化等。因此，我們有必要像設置抽象單元類一樣，設置抽象材料類，然后通過繼承的方式創(chuàng)建具體的材料類。在結構化程序設計中，大多數(shù)軟件都將材料和單元放在同一函數(shù)中處理。我們將單元和材料分開的目的是讓單元能對應不同的材料，同時一種材料能應用在不同的單元中，這種多對多的關系增加了單元和材料的應用類型。(5)荷載類這里荷載類定義的是單元荷載，其成員變量用來對荷載情況進行描述，如荷載類型、大小、作用位置和方向等。單元荷載在計算時需要先轉換成等效結點荷載，在計算單元剛度矩陣并放入總剛度矩陣時一同放入方程右端的荷載向量中，而初始結點荷載在此時直接進入方程右端的荷載向量中。三維對流擴散問題程序類設計面向對象的程序設計，歸結為類的設計。根據(jù)第三章和前面介紹的面向對象的有限元分析程序基本框架，對于三維對流擴散問題，我們只須依據(jù)具體問題在程序中增加新的類，來實現(xiàn)三維對流擴散問題有限元程序的開發(fā)。29對流擴散問題有限元分析方法類對流擴散問題有限元分析方法類是有限元方法類的派生類，它繼承了抽象方法類已經定義的所有成員變量和成員函數(shù)，主要用來控制整個程序的流程，實現(xiàn)三維對流擴散問題的有限元分析。classFConDiff:publicFinMethod{public:F_ConDiff();?F_ConDifR){};intConDiffAnlysis();//對流擴散分析程序通過調用函數(shù)intConDiffAnlysis(),就可以實現(xiàn)三維對流擴散問題有限元分析的一系列步驟。三維對流擴散問題有限元分析程序執(zhí)行順序：.執(zhí)行：CreateNodeDateChain()；〃讀入結點數(shù)據(jù)，創(chuàng)建結點鏈表.執(zhí)行：CreateMaterialChain()；〃讀入材料數(shù)據(jù)，創(chuàng)建材料鏈表.執(zhí)行：CreateElementChain()；〃讀入單元數(shù)據(jù)，創(chuàng)建單元鏈表.執(zhí)行：OutputMaterial()；outputCoor()；outputElement()；outputBC()；//初始數(shù)據(jù)輸出.執(zhí)行：Counter()；〃計算結點方程編號，確定總剛度矩陣和右端向量維數(shù).執(zhí)行：GenStiffMatrix()；〃計算單元剛度矩陣并放入總剛度矩陣.執(zhí)行：solver()；〃求解方程.執(zhí)行：outputDisplace()；//輸出結點位移以上1?3步是在三維對流擴散問題有限元分析類的構造函數(shù)F_ConDiff()中完成的，也就是說，對一個具體問題分析時，就已經將這個問題一些基本數(shù)據(jù)存放在這個方法類中了，這正符合面向對象思想——對象=數(shù)據(jù)+算法。4~8步則在voidStaticAnlysis()中完成。結點類結點類封裝了結點的信息，包括自由度、坐標、結點位移、約束條件等。對于三維對流擴散問題，我們要依據(jù)結點類的定義來輸入所有結點的相關信息。根據(jù)三維對流擴散方程組的特點，要確定結點序號、結點自由度、結點維數(shù)、結點約束條件、結點約束值、結點坐標和結點處的速度值。結點序號是輸入結點時的排列順序。結點自由度為lo結點維數(shù)為3O結點約束條件包括三種條件：一、結點屬于區(qū)域Q內點，用0來標識；二、結點在邊界B上，用1來標識；三、結點30在「上，用2來標識。再分別給出在口和「2上的約束值，另外還要輸入結點的三個坐標分量和結點處的速度在三個方向上的分量值。結點類是數(shù)據(jù)輸入的關鍵，同時對結點鏈表和單元鏈表的創(chuàng)建也起著非常重要的作用。結點信息的輸入通過成員函數(shù)NodeData(intilndjntiNodeFree,intiDimjntflage)來實現(xiàn)。NodeData::NodeData(intilndjntiNodeFree,intiDim,intflage