兩相流大作業(yè)

1、水水蒸汽兩相相變界面的數(shù)值模擬 兩相流動與熱物理大作業(yè)姓名 張蛟龍_學號 201328013524021_班級 物理308_ 指導(dǎo)教師 劉捷_ 完成時間 _2014.5.8_ 水水蒸汽兩相相變界面的數(shù)值模擬報告一文獻綜述 作為化石資源的替代產(chǎn)品，核能的高效，清潔一直備受青睞，然而光環(huán)之下，核廢料的處理不禁讓人黯然神傷。強致命性輻射，動輒千年的半衰期，慣用的辦法只能是深埋，等待下一代的聰明才智。與此同時，核廢料的利用和加速衰減一直是核能大國們的研究重點。歐洲的ADS系統(tǒng)第六代散裂靶模型計劃的目標就是要驗證高水平的核廢料轉(zhuǎn)換的可行性。散裂靶作為連接加速器和核廢料的裝置需要工作在高輻射和高熱流密度的

2、條件下，因此散裂靶的設(shè)計是ADS系統(tǒng)研制最有挑戰(zhàn)的部分。由加速器產(chǎn)生的高能質(zhì)子流轟擊靶核產(chǎn)生中子作為外源中子驅(qū)動和維持次臨界堆的運行。散裂靶在極小的空間內(nèi)需承受極大的熱負荷，質(zhì)子束通道與靶核的自由面相鄰更加劇了設(shè)計難度。受材料限制，流體的溫度不能超過550度，因此必須保證流體維持在一定的流量。但同時又要考慮高流速帶來的飛濺和回流造成的局部溫度過高。這一裝置在水作為散裂靶的實驗中獲得了成功。2 問題描述2.1.模型及尺寸 圖1、歐洲液態(tài)金屬散裂靶V0.10示意圖1如圖1所示的歐洲加速器驅(qū)動次臨界堆（ADS）之無窗散裂靶示意圖，液態(tài)鉛鉍合金從上方管間流下并匯合，形成兩相界面，質(zhì)子束由中間的真空管進

3、入打在自由面上。此次模擬用的是水，詳細物理背景見文獻1。2.2. 控制方程 連續(xù)性方程 動量方程 能量方程 3 Openfoam求解 有關(guān)Openfoam的下載和安裝在老師給的安裝指導(dǎo)的推薦網(wǎng)站上有詳細的操作，在此就不贅述。網(wǎng)址為：/download/ubuntu.php。3.1. OpenFoam求解簡述 Openfoam是一款基于linex的開源可編程軟件，其求解過程的關(guān)鍵是三個文件夾的設(shè)置，即0，constant和system。0文件夾里存放的是初始條件和邊界條件設(shè)置文件；constant文件夾里存放的是網(wǎng)格文件，物性參數(shù)和求解器模型；syst

4、em文件夾里存放的是求解過程控制，差分格式和代數(shù)方程求解器設(shè)置文件。以下就三個文件的設(shè)置展開簡述初始條件、邊界條件、物性參數(shù)，網(wǎng)格個數(shù)、疏密設(shè)置差分格式、界面捕獲算法、氣蝕模型等的選擇和設(shè)置。3.2. 0文件夾 包含有5個文件，分別為alph-water，p_rgh，U，epsilon，k，詳細設(shè)置見附錄1，這里只著重強調(diào)在大作業(yè)完成過程中幾個曾經(jīng)連續(xù)考慮的點。 首先是參數(shù)的量綱設(shè)置。在Openfoam文件中常會見到這樣一行代碼：dimensions 0 0 0 0 0，這便是量綱，單位順序依次是 質(zhì)量，長度，時間，溫度，物質(zhì)的量，電流，光強。 其次是邊界條件和初始條件的設(shè)置。在alph-wa

5、ter中，alpha代表水所占比例，參照userguide，1時表示全部為液相，0時表示全部為氣相。初始內(nèi)部場的設(shè)置均為1，即起始時刻，散裂靶內(nèi)部充滿水。水入口是邊界類型為“定值”，即fixedValue；上下兩個出口的邊界類型為“進出口”，即inletOutlet。在p_rgh中初始內(nèi)部場設(shè)置值為5330，為壓力較大的出口的壓力值，查閱資料推薦內(nèi)部場使用大值以減小汽蝕。質(zhì)子束通道理論上應(yīng)為真空，但考慮到自由面上水的汽化，選定了水的飽和壓力2330以減小汽化，保證系統(tǒng)正常工作。文獻顯示上下兩個出口壓力差選定在2000到3000之間，此次試驗選用差值為3000，即下面出口值設(shè)定為5330，也就是

6、前面提到的出場值。在U中，因為速度時矢量，在設(shè)定時需要根據(jù)網(wǎng)格的情況處理。初始內(nèi)部場設(shè)定為uniform（0 0 0），即均勻的靜止態(tài)。進水口的類型為fixedValue，且有一個1m/s的向下流的速度。兩個出口的類型為pressureInletOutletVelocity，即由壓力決定速度的邊界條件。在epsilon和k中，參照userguide中給出的公式：；， ，其中D為等效直徑，分別為散裂靶上口總的直徑和質(zhì)子入口的直徑，計算得k=0.0025，=0.00327。參照蘇軍偉的書，水入口類型為fixedValue(定值)，兩個出水口的類型為inletOutlet，而墻的類型為epsilon

7、WallFunction和kqRWallFunction。3.3. constant文件夾 首先介紹網(wǎng)格繪制。網(wǎng)格繪制有兩種方法，在大作業(yè)探索過程中都進行了嘗試。第一，關(guān)于Gambit直接繪制三維網(wǎng)格然后導(dǎo)進Openfoam進行計算。網(wǎng)格數(shù)大概為3萬多，導(dǎo)入的操作過程如下：在Openfoam中建立一個文件夾，注意此文件夾必須包括system這個文件夾，因為導(dǎo)入過程需要system中所含的controlDict文件，否則將無法導(dǎo)入，可以先“借用”其中任意一個例子的system，后期計算時再根據(jù)實際情況更改。將網(wǎng)格文件復(fù)制到該文件夾下，在控制終端輸入命令： fluentMeshToFoam ni.

8、msh -scale 0.001其中fluentMeshToFoam為導(dǎo)入命令，ni.msh為網(wǎng)格名稱，-scale 0.001的作用是把Gambit中默認為毫米的單位改成Openfoam中默認的以米為單位。單位轉(zhuǎn)換的操作對初學者尤為要注意，開始的很多次嘗試失敗的原因就是忽略了兩個軟件之間單位的不相容性。導(dǎo)入完成后會生成Polymesh的文件夾，所包含的正是全部的網(wǎng)格信息。第二，關(guān)于Gambit繪制二維網(wǎng)格后導(dǎo)入Openfoam進行旋轉(zhuǎn)。注意在Gambit繪制的二維網(wǎng)格在保存時一定要保存成二維的格式。與三位的導(dǎo)入相同，導(dǎo)入命令是： fluentMeshToFoam ni.msh -scale

9、0.001同樣生成constant文件夾然后是旋轉(zhuǎn)成五度的楔形的操作，命令是： makeAxialMesh -axis ax -wedge frontAndBackPlanes其中makeAxialMesh是生成軸網(wǎng)格的命令，-axis ax表示制定的旋轉(zhuǎn)軸是邊界名稱型為ax的邊，-wedge frontAndBackPlanes表示繞軸的旋轉(zhuǎn)面是名稱為frontAndBackPlanes的面。當然完成這一操作的前提是系統(tǒng)裝有makeAxialMesh的軟件?？梢栽诳刂平K端直接輸入makeAxialMesh這個命令，然后按提示進行下載和安裝就可以使用。命令完成后終端會提示進行網(wǎng)格清理，命令為：

10、 collapseEdges在執(zhí)行這一命令前，需要在system文件夾中添加一個collapseDict;文件規(guī)定網(wǎng)格清理的條件。具體設(shè)置見附錄。之后可以檢查一下網(wǎng)格，命令為： checkMesh這三個操作完成會生成一個新的包含網(wǎng)格信息的文件夾，它就可以作為新的0文件夾供計算使用。 然后是其余的幾個文件。包括transportProperties，turbulenceProperties，RASProperties，g。在物性參數(shù)設(shè)置transportProperties中，水的飽和壓力設(shè)定為2370Pa，表面張力sigma為0.07。流動選擇湍流模型中的不可壓縮RAS模型，選用k-兩方程模型

11、進行求解。在重力場設(shè)置時，因為也是矢量，所以要根據(jù)導(dǎo)進去的模型的坐標關(guān)系進行適當設(shè)定。具體設(shè)置見附錄。3.4. system文件夾 system文件夾中一共包含了三個文件： controlDict，fvSchemes和fvSolution分別在求解過程控制，離散格式，求解器選擇中發(fā)揮作用。首先看一下controlDict，文件中選定的是interPhaseChangeFoam，即帶有自由表面流求解器。參照數(shù)值傳熱資料的設(shè)定Co=0.5，為了保證物理過程的清晰和礙于Co數(shù)的限制，時間步長按老師要求取作1e-5。參照前期同學的經(jīng)驗，采用adjustabeRunTime，即可調(diào)步長，自動調(diào)節(jié)最后一次

12、的時間步長，以便準確輸出。開啟adjsutTimeStep，由于最大Co數(shù)的限制，當程序計算值大于最大Co數(shù)時，能夠根據(jù)最大Co數(shù)自動調(diào)節(jié)減小步長；同時開啟runtimeModifiable,在solver允許調(diào)節(jié)步長時，可根據(jù)最大Co數(shù)自動調(diào)節(jié)加大時間步長。 下面看數(shù)值格式，也就是各個量的離散形式。時間(ddtSchemes)的離散采用一階隱式Euler形式。梯度格式（gradSchemes）采用有限元高斯線性插值。對流項的離散（divSchemes）格式選取特別要注意，由數(shù)值傳熱學的知識可以知道，二階迎風有著優(yōu)良的遷移和守恒特性，所以速度的離散采用具有二階精度的二階迎風，其他參數(shù)則采用一階

13、迎風格式。擴散項（laplacianSchemes）采用二階高斯守恒格式。表面插值（interpolationSchemes）采用中心差分格式。其他的參數(shù)設(shè)置見附錄。關(guān)于界面捕獲方法，此次模擬采用的是VOF方法，即Volume of fluidmethod來捕獲界面。接下來是代數(shù)方程求解器的設(shè)定。參照Openfoam當中的例子，Alph-water的計算選用PBiCG求解器，即預(yù)條件雙共軛梯度求解器，殘差（tolerance）設(shè)定為1e-8，相對殘差（relTol）設(shè)定為0，當程序運行結(jié)果滿足兩個條件之一時停止。計算采用的是改進版的simple算法pimple，壓力場與速度場分別存放在兩套不同

14、的網(wǎng)格中，所以壓力求解器選用的是代數(shù)多重網(wǎng)格求解器（GAMG）。參照孫軍委的指導(dǎo)書，U，k，全部選用的是光滑求解器（smoothSolver）。在pimple算法設(shè)定里，壓力修正次數(shù)（nCorrectors）設(shè)為2，不進行非正交修正（nNonOrthogonalCorrectors）。其他設(shè)置見附錄。4 計算流程和結(jié)果4.1. 運行界面以下為在1.3826s時的運行程序界面：DILUPBiCG: Solving for alpha.water, Initial residual = 1.01747e-06, Final residual = 9.65653e-09, No Iterations

15、 1Phase-1 volume fraction = 0.581635 Min(alpha1) = 1.59102e-05 Max(alpha1) = 1MULES: Correcting alpha.waterLiquid phase volume fraction = 0.581612 Min(alpha1) = 1.59102e-05 Max(alpha1) = 1GAMGPCG: Solving for p_rgh, Initial residual = 0.000606079, Final residual = 5.81781e-09, No Iterations 8GAMGPCG

16、: Solving for p_rgh, Initial residual = 8.81889e-06, Final residual = 3.87137e-09, No Iterations 5smoothSolver: Solving for epsilon, Initial residual = 0.000916708, Final residual = 9.80608e-07, No Iterations 33smoothSolver: Solving for k, Initial residual = 0.000708262, Final residual = 9.89131e-07

17、, No Iterations 122ExecutionTime = 4.41 s ClockTime = 4 sCourant Number mean: 0.12362 max: 0.370546deltaT = 0.0002Time = 1.3826 我們可以看到運行程序里依次是對各個變量的求解，殘差計算，迭代次數(shù)以及Co數(shù)和時間步長。另外，我們可以發(fā)現(xiàn)程序中每次運算都會計算Liquid phase volume fraction 等，這便是界面捕獲VOF方法的體現(xiàn)。4.2. 運行結(jié)果及物理分析 以在Openfoam中旋轉(zhuǎn)生成的網(wǎng)格為例，流動在1.5s前已經(jīng)達到穩(wěn)定。 以下為alpha-w

18、ater的瞬態(tài)和穩(wěn)態(tài)圖 由圖，我們可以看出，模擬開始后首先在散裂靶的水進口末端和突擴部分的頂端由于壓力小于水的飽和壓力而產(chǎn)生了汽化現(xiàn)象。穩(wěn)態(tài)時在質(zhì)子束通道和突擴的上半部分形成了穩(wěn)定的近似真空區(qū)。以下為p_rgh的瞬態(tài)和穩(wěn)態(tài)圖上面第一個圖是流動起始時壓力的分布，可以看到，質(zhì)子束口壓力隨著水面的下降穩(wěn)定在2330Pa的設(shè)定值，水入口壓力在流場中是最高的。由第二張圖可以看出，出現(xiàn)汽化的部位的壓力為水的飽和壓力。最后的穩(wěn)態(tài)圖與alpha-water的穩(wěn)態(tài)圖幾乎一致。以下為U的瞬態(tài)和穩(wěn)態(tài)圖 由圖可以看出，在起始時刻，由于出口壓力較低，在重力作用下，突擴管道內(nèi)水的流速急劇增大，流速增大導(dǎo)致壓力減小，局部出

19、現(xiàn)汽化區(qū)，壓力下降并穩(wěn)定在飽和壓力，液體又再次流回汽化區(qū)，如此往復(fù)作用，最終在重力和壓力的聯(lián)合作用下行程了穩(wěn)態(tài)。理論上過程是這樣的，但在模擬的過程中，一直沒有得到理想的穩(wěn)態(tài)速度分布圖，還需要繼續(xù)修改一些參數(shù)。五體會啟示 老師依然是那個風格老師，作業(yè)卻是更有挑戰(zhàn)性的作業(yè)。這學期一直在學數(shù)值傳熱加上一個有關(guān)fluent的談?wù)撜n，所以對fluent已經(jīng)不是那么陌生。不過，讓人沒想到的是懷柔這一片科研凈土竟給了我們機會接觸一下開源的Openfoam。沒聽過倒不打緊，關(guān)鍵是得在linux上運行。大作業(yè)開始的前期日子，我們一直在積極探索如何安裝虛擬機和ubuntu。計算機本來內(nèi)存就已經(jīng)很小，運行虛擬機簡直

20、有些讓人抓狂，再加上老師說得巨大的計算量，我們開始考慮放棄虛擬機，改向直接裝并行系統(tǒng)。無意中得知班里有平時玩linux的高手，請來不到一個小時完成安裝。一下子感覺前些天的努力都白搭了，看來老師提醒我們的要學會利用資源是很有必要的。接下來便是繪制網(wǎng)格，前期的計劃是在較為熟悉的Gambit里繪制三維圖形然后導(dǎo)如Openfoam。考慮到計算量的問題，并沒有繪制完整的三維圖形，而是在原模型的基礎(chǔ)上截了一個五度的楔形用于模擬。網(wǎng)格繪制難度并不大，之后參照Openfoam的userguide，網(wǎng)格的導(dǎo)入同樣是比較容易。之后就進入到了艱難的初始條件設(shè)置和求解設(shè)置過程，由于參照Openfoam中的例題syst

21、em的設(shè)置大同小異，所以主要心思就放在了初始條件設(shè)置上，但在一系列嘗試都失敗后轉(zhuǎn)而懷疑是網(wǎng)格的原因，所以才有了重新學習在Openfoam中旋轉(zhuǎn)二維網(wǎng)格。網(wǎng)格的繪制總是比較簡單，信心滿滿以為這次差不多了，但運行結(jié)果還是老樣子，運行百步后Co數(shù)太大導(dǎo)致時間步長太小超過了計算機的浮點數(shù)，系統(tǒng)自動溢出。無奈之下，向理化所的王超同學請教。他認為可能是求解器的設(shè)置有問題，于是在依據(jù)他的建議進行了求解器設(shè)置的更改后計算才慢慢有了起色。以上便是本學期大作業(yè)的一個概述。 依然覺得設(shè)置大作業(yè)很有必要，并且是加大難度很有必要。因為事實證明，只要問題可解，我們總是可以找到解決的方法，只不過是前期做的功課的長短問題。首

22、先，經(jīng)過此次散裂靶的模擬，我們認識了一款Openfoam軟件，在與其斗爭的過程中，接觸到了很多內(nèi)部參數(shù)的設(shè)置選擇，這是在fluent等軟件計算時所看不到的。這樣便為我們深入探究每個參數(shù)和設(shè)定對計算的影響提供了很好的平臺嘗試機會。然后就是這個學習的過程，當接觸到一個完全陌生的軟件或者問題時先應(yīng)該從哪里下手，如何快速的找到自己需要的元素并加以利用；當在探索過程遇到障礙時如何有條理地進行分析，找到癥結(jié)所在。最后要感謝老師又給我們提供了一次“痛然后快樂”的機會，感謝王超，范俊輝和鄒文杰同學的無私幫助，感謝朱鵬同學分享的蘇軍偉老師的書，感謝陪我一起探索的欒奕軍和丁林超同學。參考文獻1 A.G. Clas

23、s, D. Angeli, A. Batta, etc. XT-ADS Windowless spallation target thermohydraulic design & experimental setup, Journal of Nuclear Materials, 415, 378-384, 2011.2 A. Batta, A.G. Class, Numerical investigations on Geometrical Designs of the Windowless XT ADS spallation target, Proceedings of ICAPP

24、2007, Nice, France, May 13-18, 2007.3 A. Batta, A.G. Class, Free surface modeling and simulation of the water experiment for the XT ADS spallation target, Proceedings of ICAPP 08, Anaheim, CA USA, June 8-12, 2008.4 F. Roelofs, N.B. Siccama, H. Jeanmart, K. van Tichelen, M. Dierckx, P. Schuurmans, Ap

25、plication of a controlled swirl in the XT-ADS spallation target, Proceedings of ICAPP 08, Anaheim, CA USA, June 8-12, 2008 說明：in1為質(zhì)子束通道的口，in2為水入口，ax為軸， frontAndBackPlanes是最初導(dǎo)進去的二維的旋轉(zhuǎn)面，frontAndBackPlanes_pos與frontAndBackPlanes_neg為旋轉(zhuǎn)后生成的楔形的前后兩個面。附錄一；初始條件設(shè)置alpha-waterdimensions 0 0 0 0 0;internalField

26、 uniform 1;boundaryField in2 type fixedValue; value $internalField; in1 type inletOutlet; inletValue uniform 1; value $internalField; out type inletOutlet; inletValue uniform 1; value $internalField; wall type zeroGradient; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos

27、 type wedge; frontAndBackPlanes_neg type wedge; epsiloninternalField uniform 0.00327;boundaryField wall type epsilonWallFunction; value uniform 0.00327; Cmu 0.09; kappa 0.41; E 9.8; value uniform 0.00327; out type inletOutlet; inletValue uniform 0.00327; value uniform 0.00327; in2 type fixedValue; v

28、alue uniform 0.00327; in1 type inletOutlet; inletValue uniform 0.00327; value uniform 0.00327; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; kinternalField uniform 0.0025;boundaryField wall type kqRWallFunction; value unif

29、orm 0.0025; out type inletOutlet; inletValue uniform 0.0025; value uniform 0.0025; in2 type fixedValue; value uniform 0.0025; in1 type inletOutlet; inletValue uniform 0.0025; value uniform 0.0025; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPl

30、anes_neg type wedge; p_rghinternalField uniform 5330boundaryField in2 type zeroGradient; in1 type totalPressure; p0 uniform 2330; U U; phi phi; rho rho; psi none; gamma 1; value $internalField; out type totalPressure; p0 uniform 5330; U U; phi phi; rho rho; psi none; gamma 1; value $internalField; w

31、all type zeroGradient; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; uinternalField uniform (0 0 0);boundaryField in2 type fixedValue; value uniform (0 -1 0); in1 type pressureInletOutletVelocity; phi phi; value $internalF

32、ield; out type pressureInletOutletVelocity; phi phi; value $internalField; wall type fixedValue; value uniform (0 0 0); ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; 附錄二：constant設(shè)置RASPropertiesRASModel kEpsilon;turbulence

33、on;printCoeffs on;transportPropertiesphases (water vapour);phaseChangeTwoPhaseMixture SchnerrSauer;pSat pSat 1 -1 -2 0 0 2370;sigma sigma 1 0 -2 0 0 0 0 0.07; water transportModel Newtonian; nu nu 0 2 -1 0 0 0 0 9.9838e-07; rho rho 1 -3 0 0 0 0 0 998.12;vapour transportModel Newtonian; nu nu 0 2 -1

34、0 0 0 0 5.5526e-04; rho rho 1 -3 0 0 0 0 0 0.017529;KunzCoeffs UInf UInf 0 1 -1 0 0 0 0 1.0; tInf tInf 0 0 1 0 0 0 0 0.5; Cc Cc 0 0 0 0 0 0 0 10; Cv Cv 0 0 0 0 0 0 0 10;MerkleCoeffs UInf UInf 0 1 -1 0 0 0 0 1.0; tInf tInf 0 0 1 0 0 0 0 0.5; Cc Cc 0 0 0 0 0 0 0 80; Cv Cv 0 0 0 0 0 0 0 1e-03;SchnerrSa

35、uerCoeffs n n 0 -3 0 0 0 0 0 1.0e+12; dNuc dNuc 0 1 0 0 0 0 0 2.0e-06; Cc Cc 0 0 0 0 0 0 0 1; Cv Cv 0 0 0 0 0 0 0 1; turbulenceProperties simulationType RASModel;附錄三：system設(shè)置controlDictapplication interPhaseChangeFoam;startTime 0;stopAt endTime;endTime 5;deltaT 5e-5;writeControl adjustableRunTime;wr

36、iteInterval 0.001;purgeWrite 0;writeFormat ascii;writePrecision 6;writeCompression uncompressed;timeFormat general;runTimeModifiable yes;adjustTimeStep on;maxCo 0.5;maxAlphaCo0.5;maxDeltaT1; fvSchemesddtSchemes default Euler;gradSchemes default cellLimited Gauss linear 1.0;divSchemes default none; div(phi,alpha) Gauss upwind; div(phirb,alpha) Gauss upwind ; div(