有限體積法基本原理 課件 chap06-2D輸運(yùn)方程及有限體積法

大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析優(yōu)化與CAE軟件全國(guó)重點(diǎn)實(shí)驗(yàn)室有限體積法基本原理夏廣慶,張軍軍,陳沖,

韓亞杰,鹿暢6.12D輸運(yùn)方程通用溫度輸運(yùn)的對(duì)流擴(kuò)散方程：忽略時(shí)間導(dǎo)數(shù)項(xiàng)和對(duì)流項(xiàng)，可以得到：在直角坐標(biāo)系下將梯度算子與點(diǎn)積算子展開(kāi)，可以得到：6.2控制方程離散控制方程中各項(xiàng)在有限體積內(nèi)積分：

6.2控制方程離散根據(jù)散度定理將控制方程中體積積分項(xiàng)進(jìn)行變換：?jiǎn)卧哂杏邢迶?shù)量的面（??），對(duì)于四邊形單元，??=4，而對(duì)于三角形單元，??=3。因此，對(duì)整個(gè)表面的曲面積分可以用對(duì)單元的??個(gè)面積分來(lái)代替：在二階有限體積法中，流動(dòng)量在整個(gè)面上都是線(xiàn)性變化的，因此，整個(gè)面的積分可以簡(jiǎn)化為面心的值（常數(shù)）：6.2控制方程離散——內(nèi)部單元離散格式為進(jìn)一步簡(jiǎn)化方程，參考圖6–1中的2D四邊形單元，該單元有左面（??）、右面（??）、頂面（??）和底面（??）。單位法向量總是指向單元外，??方向從左到右為正，??方向從下到上為正。這種離散方法適用于網(wǎng)格中所有單元形狀和方向與圖6–1中四邊形單元（正規(guī)網(wǎng)格）相同，為進(jìn)一步簡(jiǎn)化，接下來(lái)將對(duì)內(nèi)部單元和邊界單元的離散格式分別考慮。6.2控制方程離散——內(nèi)部單元離散格式所有擴(kuò)散項(xiàng)都將使用中心差法，使用圖6–2中的符號(hào)表示單元中心點(diǎn)之間的距離??，上式可以表示為：

按照單元中心溫度（????、????、????、????和????）對(duì)上式進(jìn)行排序，可以得到：6.2控制方程離散——內(nèi)部單元離散格式所以對(duì)于內(nèi)部單元，2D擴(kuò)散方程的有限體積離散格式及系數(shù)如下：使用求和符號(hào)對(duì)離散格式進(jìn)行表示：6.2控制方程離散——邊界單元（Dirichlet）右、上和下表面與內(nèi)部單元相連，因此，通過(guò)這些面的熱擴(kuò)散可采用與內(nèi)部單元相同的處理方法來(lái)進(jìn)行計(jì)算：本章將在左面應(yīng)用定值（Dirichlet）邊界條件。左邊界面上的溫度為??????????，到壁面的距離為到下一個(gè)單元中心點(diǎn)距離的一半。因此，通過(guò)左邊界面的熱擴(kuò)散為：將上式代入邊界單元對(duì)應(yīng)離散格式中：6.2控制方程離散——邊界單元（Dirichlet）為與內(nèi)部單元離散格式保持一致，將上式變換為：需要注意的是，使用求和符號(hào)時(shí)對(duì)于邊界面????

=0，并且????的求和不包括邊界面。6.2控制方程離散——邊界單元（Neumann）邊界施加固定通量（Neumann）時(shí)的邊界單元如圖6-4所示，該單元右、上和下表面與內(nèi)部單元相連，因此，通過(guò)這些面的熱擴(kuò)散可采用與內(nèi)部單元相同的處理方法來(lái)進(jìn)行計(jì)算：

6.2控制方程離散——邊界單元（Neumann）為與內(nèi)部單元離散格式保持一致，將上式變換為：需要注意的是，使用求和符號(hào)時(shí)對(duì)于邊界面????

=0，并且????的求和不包括邊界面。6.3控制方程各單元系數(shù)與線(xiàn)性系統(tǒng)根據(jù)前述推導(dǎo)，可以得到控制方程對(duì)應(yīng)單元各系數(shù)(Dirichlet邊界與內(nèi)部單元)：Neumann邊界與內(nèi)部單元系數(shù)表：?jiǎn)卧?lèi)型Interior0Left0Right0Bottom0Top0單元類(lèi)型Interior0Left00Right00Bottom00Top006.3控制方程各單元系數(shù)與線(xiàn)性系統(tǒng)將前述單元系數(shù)用求和符號(hào)表示，如下表：?jiǎn)卧骖?lèi)型Interior0Dirichlet0Neumann00

6.2控制方程離散——多邊界界單元對(duì)于2D情況，一些邊界單元可能具有不止一個(gè)邊界面，且不同邊界面的邊界條件也可能不同，如圖6-5所示，計(jì)算域左上角具有2個(gè)邊界面，該單元的右面與下面與內(nèi)部單元相連，因此，通過(guò)這些面的熱擴(kuò)散通量可采用與內(nèi)部單元相同的處理方法來(lái)進(jìn)行計(jì)算：：如果邊界面的左側(cè)為固定溫度??????????,??，頂部為固定溫度??????????,t，那么也可以在邊界面上使用中心差分法。6.2控制方程離散——多邊界界單元為與內(nèi)部單元離散格式保持一致，將上式變換為：使用求和符號(hào)表示：6.3Example1——帶熱源的熱傳導(dǎo)根據(jù)問(wèn)題模型，可以知道桿件中溫度滿(mǎn)足如下控制方程：6.3Example1——帶熱源的熱傳導(dǎo)Step-1：網(wǎng)格劃分Step-2：添加材料屬性Step-3：計(jì)算方程中各系數(shù)單元1010100-4071006021010100-2051005031010100-205100504100100-409100605

0101010-202100506

101010100100407

101010100100408

1001010-204100506.3Example1——帶熱源的熱傳導(dǎo)Step-3：計(jì)算方程中各系數(shù)單元90101010-2021005010

101010100100401110101010010040121001010-2041005013

010010-40510060141010010-2031005015

1010010-2031005016100010-40710060Step-3：計(jì)算方程中各系數(shù)(采用求和符號(hào)表示各系數(shù))單元面類(lèi)型Interior10100Dirichlet(L)2002000Dirichlet(B)2003000Dirichlet(R)2004000Dirichlet(T)20050006.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器構(gòu)造——求解器結(jié)構(gòu)及各文件功能簡(jiǎn)介(復(fù)習(xí))files文件主要完成兩個(gè)工作：1）指定需要編譯的文件，這里的文件是代碼文件，即.C文件，而不包括.H文件；2）指定編譯類(lèi)型及對(duì)應(yīng)類(lèi)型名稱(chēng)，如編譯成庫(kù)類(lèi)型，其關(guān)鍵字為L(zhǎng)IB=···；編譯成可執(zhí)行程序類(lèi)型，其關(guān)鍵字為EXE=···。options文件也主要完成兩個(gè)工作：1）指定編譯的頭文件所在目錄，如果所用頭文件在當(dāng)前文件夾或?yàn)闃?biāo)準(zhǔn)C++頭文件，無(wú)需指定，EXE_INC和關(guān)鍵字-I用于指定頭文件所在目錄；2）指定編譯當(dāng)前程序用到的庫(kù)，EXE_LIBS和關(guān)鍵字-L/-l用于指定用到的庫(kù)。6.4Example1——帶熱源的熱傳導(dǎo)createFields.H——初始化基本物理場(chǎng)待求溫度場(chǎng)速度場(chǎng)溫度擴(kuò)散系數(shù)熱源OpenFOAM求解器構(gòu)造OpenFOAM求解器構(gòu)造6.4Example1——帶熱源的熱傳導(dǎo)convectionAndDiffusion.C——求解器源文件網(wǎng)格讀取線(xiàn)性方程構(gòu)造方程求解及輸出算例文件路徑正確性檢測(cè)時(shí)間類(lèi)初始化網(wǎng)格類(lèi)初始化基本物理場(chǎng)初始化構(gòu)造線(xiàn)性方程：各項(xiàng)與算例問(wèn)題中對(duì)應(yīng)線(xiàn)性方程求解結(jié)果文件輸出6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器構(gòu)造編輯完成OpenFOAM各求解器后，在求解器下執(zhí)行wmake對(duì)求解器進(jìn)行編譯，得到可執(zhí)行程序。求解器源文件編譯得到求解器名稱(chēng)g++之后的指令參數(shù)依次指定了采用的C++標(biāo)準(zhǔn)版本、編譯器參數(shù)、機(jī)器位數(shù)、整數(shù)位數(shù)等信息，之后以“-I”開(kāi)頭的表示使用到的頭文件所在的目錄，以“-L”開(kāi)頭的表示指定的庫(kù)，以“-1”開(kāi)頭的表示鏈接到指定的庫(kù)，“-oapp”表示生成的程序名為app。6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM

求解器對(duì)應(yīng)算例構(gòu)造——帶熱源熱傳導(dǎo)算例構(gòu)造算例文件較為方便方法修改已有的、與問(wèn)題模型相似的算例，這里以修改cavity為例進(jìn)行介紹，對(duì)應(yīng)算例文件路徑：openfoam10/tutorials/incompressible/icoFoam/cavity6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器對(duì)應(yīng)算例構(gòu)造——帶熱源熱傳導(dǎo)算例6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器對(duì)應(yīng)算例構(gòu)造——帶熱源熱傳導(dǎo)算例6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器對(duì)應(yīng)算例構(gòu)造算例文件編輯完成后以此執(zhí)行：blockMesh、../app、paraFoam分別進(jìn)行網(wǎng)格劃分、算例求解及后處理，對(duì)應(yīng)各步驟過(guò)程分別如下：6.4Example1——帶熱源的熱傳導(dǎo)OpenFOAM求解器對(duì)應(yīng)算例構(gòu)造算例文件編輯完成后以此執(zhí)行：blockMesh、../diffusionEquationSolver、paraFoam

分別進(jìn)行網(wǎng)格劃分、算例求解及后處理，對(duì)應(yīng)各步驟過(guò)程分別如下：6.4Example1——帶熱源的熱傳導(dǎo)自編程序簡(jiǎn)介——2D有熱源熱傳導(dǎo)通過(guò)自編程序方式進(jìn)行問(wèn)題求解，程序主要包含模型參數(shù)設(shè)置、網(wǎng)格劃分、矩陣系數(shù)計(jì)算、線(xiàn)性系統(tǒng)構(gòu)造、求解及后處理。程序各部分定義參考如下：6.4Example1——帶熱源的熱傳導(dǎo)自編程序簡(jiǎn)介——2D有熱源熱傳導(dǎo)通過(guò)自編程序方式進(jìn)行問(wèn)題求解，程序主要包含模型參數(shù)設(shè)置、網(wǎng)格劃分、矩陣系數(shù)計(jì)算、線(xiàn)性系統(tǒng)構(gòu)造、求解及后處理。程序各部分定義參考如下：6.4Example1——帶熱源的熱傳導(dǎo)算例結(jié)果在進(jìn)行后處理時(shí)較為方便的繪圖方法是將結(jié)果保存在網(wǎng)格節(jié)點(diǎn)上而非網(wǎng)格單元中心，如圖6–10所示，因?yàn)檫吔缟系臏囟戎凳欠植荚诰W(wǎng)格節(jié)點(diǎn)上（圖中紅色部分），如果采用單元中心（圖中藍(lán)色部分）進(jìn)行后處理則無(wú)法得到邊界上的溫度分布。常見(jiàn)的CFD程序采用單元中心進(jìn)行計(jì)算，因此，內(nèi)部節(jié)點(diǎn)值（圖6–10中以白色顯示）需要通過(guò)單元中心點(diǎn)值（藍(lán)色顯示）之間的插值計(jì)算得到。6.5熱平衡對(duì)于2D熱平衡檢驗(yàn)時(shí)，內(nèi)部單元一般有限體積離散格式如下：根據(jù)傅里葉定律通過(guò)每個(gè)單元面的熱通量為：將上式代入對(duì)應(yīng)離散格式中可以得到：負(fù)號(hào)是為了確保熱量流動(dòng)方向與溫度梯度方向相反對(duì)于2D四邊形單元而言，熱平衡誤差表示為：6.5熱平衡計(jì)算前述算例問(wèn)題中各單元的熱平衡誤差，結(jié)果如下表所示：6.5熱平衡單元熱通量平衡表還可用于計(jì)算平板每個(gè)面（跨越所有邊界）流出的總熱通量：將每個(gè)邊界面（上表中用紅色標(biāo)出的部分）流出的熱通量相加，結(jié)果如表6–8所示：正熱通量表示熱量穿過(guò)平板從平板邊界流出，負(fù)熱通量則表示熱量通過(guò)邊界進(jìn)入平板內(nèi)。左邊界作為溫度最低的面（????

=100℃），該平面的熱通量最大。相反，上表面作為最熱面（????

=250℃），熱量通過(guò)該平面進(jìn)入平板內(nèi)，從而產(chǎn)生