【文章內(nèi)容簡介】 （325）JJJAB???? （326）JJnDT式中：△ε —— 一個時間步長內(nèi)凝固界面的運動距離/m；△rs —— 一個時間步長內(nèi)升華界面的運動距離/m。 TDMA算法計算過程中，沿徑向劃分在液相區(qū)內(nèi)的網(wǎng)格數(shù)為 75，固相區(qū)內(nèi)的網(wǎng)格數(shù)為 75，液滴內(nèi)的網(wǎng)格總數(shù)為 150，時間步長△t 為 107 s，數(shù)值求解方法采用三角對陣法（ TDMA） 。即對于網(wǎng)格離散后控制方程的統(tǒng)一形式： （ J=1~150） （327）111JJJJnnnATBCTD????顯然當(dāng) J=1 時，C J=0；而 J=150 時，B J=0，即首位兩個節(jié)點的方程中僅有兩個未知數(shù)，下面采用消元和回代的方法。消元時，從 J=2 開始，逐一把每個式子中的非零元素消去一個，使原來的三元方程化為二元方程。消元進行到最后一行時，該二元方程化為一元方程，可立即得出該未知量的解，然后逐一往回代，由各二元方程解出其它未知量。消元的目的是把（327）式化成下列形式的方程： （328）11????JJJQTP其中：13 （329）1??JJPCAB （330）1??JJQD這兩個計算 ， 的通式是遞歸的，即要計算 ， ，需要知道 ， ，最終要知JPJPQ1?JPJQ道 P1，Q 1 之值。P 1，Q 1 可以由左端點的離散方程來確定： 102DTCBA??其中： 01TC所以： ，P?1當(dāng)消元進行到最后一行時，有：，150150QT?而： ?P所以： = 150逐個往回回代，可以得出 TJ(J=149,……1)。如上所述，應(yīng)用 fortran 語言程序進行計算，首先需要對模型進行坐標(biāo)變換以解決移動界面的問題，然后對坐標(biāo)變換后的微分方程采用控制容積積分法進行網(wǎng)格離散，并用TDMA 算法進行迭代計算。最終得到某一工況下（即某一個空氣流速和液滴尺寸下） ，液滴在凝固過程中液滴溫度分布和變化情況。具體程序設(shè)計步驟為：（1） 結(jié)冰開始時刻，已知升華界面與凝固界面溫度一致 Td,s0 = Tf0。（2） 由式（210）可獲得液滴表面過冷引起的凝固界面運動速度 Vf，以及一個時間步長內(nèi)凝固界面的移動距離△ε = Vf△t。（3） 由已知空氣流速，獲得雷諾數(shù) Re、施密特數(shù) Sc，并根據(jù)公式（216）計算得到舍伍德數(shù) Sh 和傳質(zhì)系數(shù) hm。根據(jù)液滴表面溫度得到液滴表面冰層水蒸氣壓，由（27）和（215）式，得到升華界面運動速度 VS、液滴表面?zhèn)髻|(zhì)換熱的熱流密度 qm，以及一個時間步長內(nèi)升華界面的移動距離△r s = VS△t。（4） 通過坐標(biāo)轉(zhuǎn)換后的能量方程和定解條件，求解液滴內(nèi)部溫度分布，得到新的液滴表面溫度 Td,s1 和新的凝固界面溫度 Tf1。14（5） 采用新的液滴表面溫度和凝固界面溫度，重復(fù)步驟（2）~（4） ，直到滿足收斂條件： (331)11d,s,sff0., iiTT???????（6） 液滴表面溫度和凝固界面溫度均收斂后，進入下一個時間步長，重復(fù)步驟（2）~ （5） ，液滴表面半徑隨時間的變化為： (332)ss0S()trtVd???液滴內(nèi)凝固界面半徑隨時間的變化為： (333)0f()tt?? Fortran語言簡介FORTRAN，亦譯為福傳，是英文“FORmula TRANslator”的縮寫，譯為“ 公式翻譯器”，它是世界上最早出現(xiàn)的計算機高級程序設(shè)計語言，廣泛應(yīng)用于科學(xué)和工程計算領(lǐng)域。FORTRAN 語言 以其特有的功能在數(shù)值、科學(xué)和工程計算領(lǐng)域發(fā)揮著重要作用。1951 年 ， 美 國 IBM 公 司 約 翰 貝 克 斯 (John Backus)針 對 匯 編 語 言 的 缺 點 著 手 研究 開 發(fā) FORTRAN 語 言 ， 并 于 1954 年 在 紐 約 正 式 對 外 發(fā) 布 。 約 翰 貝 克 斯 提 出 的FORTRAN 語 言 稱 為 FORTRANⅠ ， FORTRANⅠ 雖 然 功 能 簡 單 ， 但 它 的 開 創(chuàng) 性 工 作 ，在 社 會 上 引 起 了 極 大 的 反 響 。 到 1957 年 第 一 個 FORTRAN 編 譯 器 在 IBM704 計 算 機上 實 現(xiàn) ， 并 首 次 成 功 運 行 了 FORTRAN 程 序 。 1958 年 ， FORTRANⅠ 進 行 了 擴 充 和 完 善 ， 引 進 了 子 函 數(shù) 等 概 念 ， 推 出 了 商 業(yè) 化的 FORTRANⅡ 版 本 。 之 后 ， FORTRAN 語 言 發(fā) 展 迅 速 ， 多 種 版 本 相 繼 在 其 它 計 算 機上 實 現(xiàn) 。Fortran 語言的主要特性有：Fortran 語言最接近數(shù)學(xué)公式的自然描述，在計算機里具有很高的執(zhí)行效率。 　　易學(xué)，語法嚴謹。 　　可以直接對矩陣和復(fù)數(shù)進行運算，這一點類似 matlab。 　　自誕生以來廣泛地應(yīng)用于數(shù)值計算領(lǐng)域，積累了大量高效而可靠的源程序。 　　很多專用的大型數(shù)值運算計算機針對 Fortran 做了優(yōu)化。 廣泛地應(yīng)用于并行計算和高性能計算領(lǐng)域。 Fortran 語 言 是 一 種 極 具 發(fā) 展 潛 力 的 語 言 ， 在 全 球 范 圍 內(nèi) 流 行 過 程 中 ， Fortran 語言 的 標(biāo) 準(zhǔn) 化 不 斷 吸 收 現(xiàn) 代 化 編 程 語 言 的 新 特 性 ， 雖 然 沒 有 現(xiàn) 在 流 行 的 的 JAVA， C++等 高 級 語 言 具 有 創(chuàng) 造 性 ， 但 是 很 多 優(yōu) 秀 的 工 程 計 算 軟 件 都 是 運 用 Fortran 語 言 編 寫 ，在 工 程 計 算 領(lǐng) 域 仍 然 占 有 重 要 地 位 。15針對本文的程序計算，下面選取幾段本文中所用到的 Fortran 程序語言對于節(jié)點的選取，程序如下： dx=1./(m1) x(1)=0. do 12 i=2,m x(i)=(i1)*dx 12 continue dr=1./(nm) do 13 i=m+1,n r(i)=(im)*dr13 continue 利用控制容積法對于模型的差分，以便于下一步的計算。程序如下：do 3000 i=1,n if() then bb(i)=1.cc(i)=0. aa(i)=1. dd(i)=0. end if if(i .ge. 2 .and. i .le. m1) then bb(i) = xkl/rhol/xcpl*dtd/ylsold/ylsold/dx amp。 +xkl/rhol/xcpl*dtd/x(i)/ylsold/ylsold amp。 +x(i)/2./ylsold*dyls cc(i) = xkl/rhol/xcpl*dtd/ylsold/ylsold/dx amp。 xkl/rhol/xcpl*dtd/x(i)/ylsold/ylsold amp。 x(i)/2./ylsold*dyls aa(i) = dx+bb(i)+cc(i) dd(i) = dx*T(i)16 end if if(i .eq. m) then bb(i) = xks/(Roldylsold)/dr cc(i) = xkl/ylsold/dx dd(i) = rhol*xclL*Viold aa(i) = bb(i)+cc(i) end if if(i .ge. m+1 .and. i .lt. n) then bb(i) = xks/rhos/xcps*dtd/(Roldylsold)/(Roldylsold)/dr amp。 +xks/rhos/xcps*dtd/((Roldylsold)*r(i)+ylsold)/(Roldylsold) amp。 +(1r(i))/2./(Roldylsold)*dyls+r(i)/2./(Roldylsold)*dysv cc(i) = xks/rhos/xcps*dtd/(Roldylsold)/(Roldylsold)/dr amp。 xks/rhos/xcps*dtd/((Roldylsold)*r(i)+ylsold)/(Roldylsold) amp。 (1r(i))/2./(Roldylsold)*dylsr(i)/2./(Roldylsold)*dysv aa(i) = dr+bb(i)+cc(i) dd(i) = dr*T(i) end if if(i .eq. n) then bb(i) = 0. cc(i) = xks/ysvold/dr aa(i) = rhos*VR*xcpv/BTnew amp。 +xks/dr/ysvold ! aa(i) = ++xks/dr/ysvold 17 dd(i)=rhos*VR*xcpv*Taold/BTnew amp。 rhos*xclS*VR! dd(i) = *Ta+*xclS*VR end if6 format(i6,)3000 continue 計算結(jié)果應(yīng)用 Fortran 程序運算我們最終得到了在一定的空氣流速和液滴尺寸下，液滴在冷空氣流中的凝固過程，液滴在不同時刻的中心溫度和表面溫度。例如當(dāng)空氣的流速為 ,液滴半徑為 時，經(jīng)過 fortran 語言的運算得到的計算結(jié)果為：表 31 fortran 語言的計算結(jié)果通過上述計算結(jié)果，利用 軟件將結(jié)果繪制到坐標(biāo)軸上，可以得到液滴中心溫度和表面溫度隨時間變化的圖形，形象給出了，液滴溫度隨著結(jié)冰過程進行的變化趨勢，通過這些坐標(biāo)圖，方便比較獲得影響規(guī)律。同時我們可以通過修改程序中空氣流速和液滴尺寸的賦值，得到不同工況下液滴在凝固時中心溫度和表面溫度的變化，方便對18比，獲得影響冷空氣流中液滴凝固的影響因素。194 計算結(jié)果分析本章內(nèi)容應(yīng)用上一章網(wǎng)格離散結(jié)果，并應(yīng)用 Fortran 語言程序，對冷空氣中液滴的結(jié)冰過程進行了計算結(jié)果分析。本文模型較 Hindmarsh 文獻模型的優(yōu)點是：引入液體過冷凝固的動力條件，并結(jié)合表面冰層升華的傳熱傳質(zhì)過程，能準(zhǔn)確預(yù)測液滴的完全凝固時間；同時對液滴內(nèi)部的溫度分布有更全面的掌握。本章計算結(jié)果分析內(nèi)容主要包括：空氣流速對液滴溫度變化的影響、液滴尺寸對其溫度變化的影響、以及液滴中心溫度與表面溫度的比較。應(yīng)用 Fortran 程序，計算結(jié)果顯示液滴結(jié)冰過程的溫度變化分為三個階段。 （1）結(jié)冰開始階段，由于凝固潛熱的釋放，液滴溫度迅速上升；（2）凝固階段，隨著凝固過程的進行，凝固界面由表面向中心運動，直至液滴完全凝固，該階段液滴溫度維持約0oC；（ 3）升華階段，液滴完全結(jié)冰后，隨