基于斐波那契數(shù)積分計算非等徑顆粒視角系數(shù)

包 濤， 高若晗， 譚援強*

(1.華僑大學(xué) 制造工程研究院，廈門 361021；2.脆性材料產(chǎn)品智能制造技術(shù)國家地方聯(lián)合工程研究中心，廈門 361021)

1 引 言

粉末床在化學(xué)工程、核工業(yè)、粉末冶金和激光增材制造等許多科學(xué)和工業(yè)領(lǐng)域中起著重要作用。粉末床在實際工況中涉及許多物理現(xiàn)象，包括顆粒材料的流動和壓實、顆粒之間的能量傳遞以及顆粒的粘結(jié)和偏析機制[1-3]。粉末床傳熱對零件成形質(zhì)量具有重要影響，引起了眾多學(xué)者的關(guān)注[4-6]。熱傳導(dǎo)、熱對流和熱輻射都有助于熱傳遞，尤其熱輻射在高溫下具有重要意義。如處于1100 ℃環(huán)境的碳化硅粉體體系內(nèi)，以熱輻射方式進(jìn)行的熱通量占據(jù)熱傳遞總量的35%；如果將材料換為玻璃顆粒，其占比將達(dá)到85%左右[7]。在顆粒燒結(jié)和熔化時，由于是局部輸入能量，粉末床中的熱輻射更加重要。

粉末床中填充顆粒的多尺度使得監(jiān)測變得極為困難，因此，很少有人在顆粒尺度進(jìn)行熱傳導(dǎo)實驗研究。通過數(shù)值模擬的方法可以更好地理解粉末床的傳熱機理[8,10]。目前，已經(jīng)提出了基于有限元和離散元等許多顆粒間的傳熱模型，還針對粉末床的有效導(dǎo)熱系數(shù)提出了幾種經(jīng)驗公式進(jìn)行預(yù)測[11,13]。然而，在大規(guī)?；騽討B(tài)顆粒系統(tǒng)中，顆粒之間的熱輻射不能忽略。但現(xiàn)有熱輻射模型在計算成本與準(zhǔn)確性這兩個方面不能同時達(dá)到數(shù)值模擬的要求。

粉末床熱輻射研究常用的模擬方法可分為兩類。第一類是基于連續(xù)體的方法，將粉末床視為多孔介質(zhì)，使用微分方程對邊界條件進(jìn)行描述。有研究表明，由于熱輻射強度的測量和計算僅在接近法線方向上一致，基于連續(xù)體不能很好地對含有大顆粒的非等溫粉末輻射特性進(jìn)行模擬[14]。第二類是基于離散的方法，將熱輻射視為顆粒表面之間發(fā)生的局部效應(yīng)，這樣，分析粉末床中熱輻射的關(guān)鍵就轉(zhuǎn)化為計算兩個顆粒之間的熱輻射通量。

采用Voronoi圖方法研究和分析粉末床中的熱輻射，該方法考慮了包括熱輻射在內(nèi)的三種傳熱機制[15-17]。但該方法的問題是無法對顆粒之間存在障礙物時的視角系數(shù)進(jìn)行計算。使用射線追蹤和蒙特卡羅等隨機的方法和其他技術(shù)也可以對粉末床中的輻射特性進(jìn)行計算[10,18-20]，但計算結(jié)果的精度在很大程度上取決于所使用的射線數(shù)或采樣點的分布，因此難以獲得高精度的結(jié)果，而且要耗費很大的計算成本。

熱輻射數(shù)值模擬的主要難點在于物體間視角系數(shù)的計算。卜昌盛等[21]對顆粒間的傳熱進(jìn)行了詳細(xì)的分析，并建立了一個考慮顆粒內(nèi)部導(dǎo)熱、顆粒粗糙表面?zhèn)鳠岷皖w粒表面氣膜與接觸間氣膜導(dǎo)熱等多種因素的傳熱數(shù)值模型。陳宇杰[22]針對顆粒間熱輻射受到阻礙這一特征，提出了考慮阻礙存在的顆粒間傳熱輻射的影響。該方法將阻礙物之間的復(fù)雜輻射情況進(jìn)行簡化，逐一計算了顆粒間的輻射熱通量。Feng等[23,24]提出了一種高效和高精度的計算方法，計算由相同大小球體組成的隨機粉末床中球體之間的視角系數(shù)。該方法的局限性主要在于沒有考慮粉末床中的顆粒粒徑的差異。大多數(shù)工業(yè)應(yīng)用中涉及的顆粒粒徑不是單一的，通常遵循高斯分布、雙峰分布和多元分布。

本文是在Feng等[23]工作的基礎(chǔ)上，提出一種新的數(shù)值程序，用于計算粉末床中非等徑顆粒之間的視角系數(shù)。首先采用斐波那契數(shù)列對非等徑顆粒表面進(jìn)行離散，并且在z方向上進(jìn)行非均勻變換以提高計算精度，最后對計算精度和效率進(jìn)行評估，并與現(xiàn)有文獻(xiàn)的結(jié)果進(jìn)行比對。

2 非等徑顆粒間無障礙視角系數(shù)計算

2.1 有限曲面間視角系數(shù)的計算

圖1 兩個有限曲面間視角系數(shù)示意圖

表面A1對表面A2的視角系數(shù)F1 - 2表示為

(1)

同理，也可以依據(jù)視角系數(shù)的定義，推導(dǎo)出表面A2對表面A1的視角系數(shù)為

F2- 1=(A1/A2)F1- 2

(2)

2.2 無阻礙顆粒非等徑顆粒間視角系數(shù)的計算

兩表面間視角系數(shù)的計算公式是基于笛卡爾坐標(biāo)系下的定積分函數(shù)。除一些簡單的幾何情況外，求解視角系數(shù)函數(shù)的解析解是極具挑戰(zhàn)性的。

因此，本文使用了一個新的兩顆粒間視角系數(shù)的計算公式，不僅減少了計算時間，還提高了計算精確度，與以往文獻(xiàn)比較驗證了該方法的可實現(xiàn)性和高效性。

假設(shè)一對半徑分別為r1和r2的球形顆粒(r1≥r2)，兩顆粒間距可表示為g≥0。兩顆粒之間的視角系數(shù)可以簡化為一個半徑為1和一個半徑為r2/r1兩顆粒間的視角系數(shù)，此時顆粒間距轉(zhuǎn)換為g/r1。假設(shè)p1=r2/r1,p2=g/r1，可將討論范圍轉(zhuǎn)換為一對半徑分別為1和p1，以及其間距為p2的一對顆粒之間的視角系數(shù)。在全局笛卡爾坐標(biāo)下，顆粒S1的球心與坐標(biāo)原點重合，顆粒S2的球心位于坐標(biāo)點(0,0,1+p1+p2)。

為了計算顆粒S1對顆粒S2的視角系數(shù)，兩顆粒分別沿z軸方向劃分為三個區(qū)域，如圖2所示。

圖2 兩顆粒中不同區(qū)域的劃分

顆粒S1對顆粒S2的視角系數(shù)可表示為[23]

F=FS - 1+FS - 2

(3)

其中

(4)

(5)

(6)

通過求解定積分可以得到兩顆粒間的視角系數(shù)的精確解，并且計算結(jié)果由兩個自變量決定，即

F=F(p1,p2)

(7)

2.3 非等徑顆粒的斐波那契積分

使用斐波那契數(shù)列原理對顆粒表面進(jìn)行離散。此方法是將顆粒放置于圓柱坐標(biāo)系中，在[0,2π]×[-r,r]區(qū)域內(nèi)將顆粒表面進(jìn)行離散。因此，兩顆粒間視角系數(shù)可以表示為

g(xi,xj)FN(p1,p2)

(8)

式中NF為斐波那契數(shù)，局部圓柱坐標(biāo)系下的第NF +1個離散點的坐標(biāo)和加權(quán)分別表示為

wi=4πr2/NF

(i=0,…，NF) (9)

(10)

(11)

式中

(12)

2.4 斐波那契積分的非均勻變換

在對兩顆粒間視角系數(shù)計算時，除采用斐波那契數(shù)列原理對顆粒表面進(jìn)行離散外，還引入非均勻變換以提高計算精度。

定義x∈[-a,a]的函數(shù)f(x)，其定積分形式為

(13)

引入新的變量x′，x′和x的關(guān)系可以定義為

x=τ(x′),x′∈[-a,a]

(14)

則式(13)可轉(zhuǎn)換為

(15)

上述非均勻變換方法應(yīng)用于z方向積分，則修正后的變換方程τ(z)可表示為

τ(z)=z+rsin(πz/r)/π

(16)

式中r為球形顆粒的半徑，對τ(z)求導(dǎo)可得

τ′(z)=1+cos(πz)/r

(17)

因此，計算視角系數(shù)的張量形式可寫為

(18)

2.5 斐波那契積分方案的計算與驗證

斐波那契數(shù)列積分方案的計算精度取決于使用的斐波那契數(shù)列項數(shù)的取值。本算例NF,p1和p2的取值分別為NF={13,21,34,55,89,144,233,377,610},p1={0.1,0.3,0.7,1},p2={0,0.01,0.02,0.05,0.1,0.5,1,2,4,6 8}。將計算結(jié)果與式(1)計算的精確值進(jìn)行比較。在一定參數(shù)p1,p2和FN條件下，斐波那契數(shù)列積分方案估算的兩顆粒間的視角系數(shù)FN(p1,p2)與精確值F(p1,p2)之間的相對誤差如圖3所示。

圖3 不同p1和p2下，斐波那契數(shù)列積分獲得兩顆粒間視角系數(shù)的相對誤差

可以看出，當(dāng)FN取為55時，已經(jīng)可以獲得較低的相對誤差值(最大5%)。

將p1=1作為特例，把計算結(jié)果和文獻(xiàn)進(jìn)行對比。如圖4所示，明顯看出當(dāng)FN=55時，兩種方法計算結(jié)果非常吻合。

圖4 斐波那契數(shù)列積分法計算獲得視角系數(shù)與文獻(xiàn)對比[22,23]

3 非等徑顆粒間存在阻礙顆粒視角系數(shù)的計算

實際顆粒體系中，兩顆粒間可能會存在多個阻礙顆粒。本節(jié)考慮了兩顆粒間存在第三個顆粒阻礙其視角情況下視角系數(shù)的計算。將粉末床中的所有顆粒視為兩目標(biāo)顆粒的障礙將耗費大量的計算成本,因此有必要減少搜索區(qū)域以更有效地檢測障礙物。球心位于陰影搜索區(qū)域的顆粒有可能成為兩個顆粒S1和S2的視角系數(shù)的阻礙顆粒，其中rmax表示顆粒體系中最大的顆粒半徑，如圖5所示。

圖5 兩顆粒間搜索中間阻礙顆粒的搜索范圍

考慮兩個顆粒S1和S2間存在一個阻礙顆粒S3情況下其視角系數(shù)的估算,只需要添加判斷兩個樣點連成的光線是否與阻礙顆粒表面相交即可。

3.1 光線阻礙檢查

在圖6中，假設(shè)光線從顆粒S1表面上樣點x1指向顆粒S2表面上樣點x2，中間阻礙顆粒的球心坐標(biāo)為O3，則阻礙顆粒球心距離光線之間的距離為

圖6 兩顆粒間視角系數(shù)(中間存在一個阻礙顆粒)

dc=‖e×r‖/‖e‖

(19)

若dc

由此，提出視角阻礙函數(shù)vb，添加第三顆粒S3后，兩顆粒間視角系數(shù)FN(p1,p2,Ω)表示為

v(xi,xj)[g(xi,xj)+

(20)

式中Ω為阻礙顆粒，且視角阻礙函數(shù)vb可表示為

(21)

3.2 計算和優(yōu)化

為驗證斐波那契數(shù)列積分法在兩顆粒間存在第三顆粒情況下估算的適用性，本文計算了當(dāng)NF={34,55,89,144,233,377,610,987}情況下兩顆粒間的視角系數(shù)。為了對估算結(jié)果進(jìn)行系統(tǒng)化對比，設(shè)定第三顆粒rc=p1，并假設(shè)第三顆粒始終與其他兩顆粒相切接觸。為方便討論，規(guī)定F(p1,p2,Ω)表示受第三顆粒阻礙后獲得的視角系數(shù)精確值，F(xiàn)(p1,p2)表示中間無阻礙顆粒的視角系數(shù)精確值，F(xiàn)(Ω)表示受第三顆粒阻礙而減少的視角系數(shù)精確值，則

F(p1,p2,Ω)=F(p1,p2)-F(Ω)

(22)

在斐波那契數(shù)列積分估算中，受第三顆粒阻礙后獲得的視角系數(shù)估算值、中間無阻礙顆粒的視角系數(shù)估算值、受第三顆粒阻礙而減少的視角系數(shù)估算值分別表示為FN(p1,p2,Ω),F(p1,p2)和FN(Ω)，則

FN(p1,p2,Ω)=FN(p1,p2)-FN(Ω)

(23)

當(dāng)FN值取足夠大時，估算值誤差極小。因此本文假設(shè)當(dāng)FN=987時獲得的兩顆粒間存在阻礙顆粒情況下獲得的估算值FN(p1,p2,Ω)即為精確值F(p1,p2,Ω)。

圖7顯示了當(dāng)參數(shù)NF={34,55,89,144,233,377,610}、p1={1,0.7,0.3,0.1}且p2={0,0.01,0.02，0.05,0.1,0.5,1,2}情況下，斐波那契數(shù)列積分法估算的兩顆粒間視角系數(shù)FN(p1,p2,Ω)與準(zhǔn)確值FN=987(p1,p2,Ω)之間的相對誤差?？梢钥闯?，相對誤差是不穩(wěn)定的，但是當(dāng)FN=55 時，其相對誤差始終小于6%。

圖7 當(dāng)rc=p1且p1=1，0.7，0.3，0.1情況下，斐波那契數(shù)列積分法估算兩顆粒間視角系數(shù)相對誤差隨p2的變化

為提高顆粒間視角系數(shù)計算的精度，需要對其進(jìn)行優(yōu)化。假定斐波那契數(shù)列積分法估算兩顆粒間視角系數(shù)在中間有無阻礙顆粒兩種情況下的估算精度相同，即

(24)

則經(jīng)過優(yōu)化后的視角系數(shù)可表示為F′(p1,p2,Ω)，即

(25)

鑒于當(dāng)NF值或p2值較小時，利用斐波那契數(shù)列積分法估算兩顆粒間視角系數(shù)精度低，對第三顆粒直徑rc=p1情況進(jìn)行優(yōu)化。求NF={34,55}，p1={1,0.7,0.3,0.1}，p2={0,0.01,0.02,0.05,0.1}情況下的視角系數(shù)的優(yōu)化值F′(p1,p2,Ω)。在NF={34,55}且p2={0,0.01,0.1}情況下，經(jīng)過優(yōu)化后的顆粒間視角系數(shù)的最大相對誤差分別由 18.11% 下降到 7.41%，5.71% 下降到 2.71%，-38.11% 下降到20.9%，12.55%下降到10.2%，如圖8所示。由此可以發(fā)現(xiàn)，在NF較小情況下，經(jīng)過優(yōu)化后的視角系數(shù)估算值，即F′(r1,r2,g,Ω)可以很大程度上提高計算精度。

圖8 當(dāng)rc=p1時，改進(jìn)后計算方法與優(yōu)化前計算方法計算精度對比

將以上提出的斐波那契數(shù)列積分法優(yōu)化后兩顆粒間視角系數(shù)與文獻(xiàn)的計算結(jié)果在圖9中進(jìn)行對比，可以看出，當(dāng)NF=55且p1=rc情況下，兩種方法獲得的估算基本一致。

圖9 優(yōu)化后斐波那契數(shù)列積分法的兩顆粒間視角系數(shù)與文獻(xiàn)計算結(jié)果對比[23,24]

4 非均勻顆粒視角系數(shù)的計算與評價

斐波那契積分方案可以進(jìn)行光線跟蹤并且對位置分布進(jìn)行優(yōu)化。為了保持合理的計算效率，需要確定球體與目標(biāo)球體之間的距離，以便在填充床中對其進(jìn)行視角系數(shù)的計算。

為了保持合理的數(shù)值精度，令NF={34,55,89,144,233,377,610}并用于后續(xù)計算。本文將通過幾個算例對非均勻顆粒間的視角系數(shù)進(jìn)行計算，并對計算的準(zhǔn)確性和效率進(jìn)行評估。

4.1 精度評估

本文利用斐波那契數(shù)列積分法計算了堆積密度分別為0.74，0.567，0.474，0.428和0.374的顆粒體系。其中第一種為標(biāo)準(zhǔn)面心立方結(jié)構(gòu)(FCC)、其余四種為修改后的面心立方結(jié)構(gòu)，即顆粒半徑r∈[rlow,rhi]的均勻分布。

選擇每個部件中心的球體作為目標(biāo)球體，球體的半徑設(shè)置為1，計算該球體與其他球體的視角系數(shù)。每個組件產(chǎn)生的空間足夠大，因此目標(biāo)球體發(fā)出的輻射將由組件中的一些球體完全接收。顯然，隨著距離目標(biāo)球體的總和距離的增加，所有視角系數(shù)的總和將趨向于1，這是總和的精確值，如圖10所示。

圖10 距離目標(biāo)顆粒不同區(qū)域半徑內(nèi)，目標(biāo)顆粒對所有顆粒視角系數(shù)的積分值

此外，接近1的速率隨顆粒分布而變化。隨著堆積密度的降低，視角系數(shù)總和趨向1的速率減慢。通常，有效區(qū)域半徑值越大，視角系數(shù)的精度越高，但涉及的計算成本越高?？梢钥闯觯?dāng)D大約為6時，顆粒視角系數(shù)的積分可以達(dá)到95%。結(jié)果表明，在不同球數(shù)的填充床中，用該方法計算的視角系數(shù)在6D內(nèi)可達(dá)到5%的相對誤差。

4.2 效率評估

為了評估提出的顆粒組合方法的整體效率，考慮了六個隨機堆積的組件，在0.8～1.2之間均勻分布。每個組件的顆粒數(shù)N分別為10，50，100，500和1000。為了進(jìn)行適當(dāng)?shù)挠嬎阈时容^，本文將每個球體的有效域D的半徑設(shè)置為6。

表1列出了兩種方法計算視角系數(shù)所需的時間。將本文方法的計算時間與蒙特卡洛法的計算時間進(jìn)行比較，計算效率提升了大約30%。

表1 含有N個顆粒的顆粒系統(tǒng)計算所有視角系數(shù)所耗費的時間

5 結(jié) 論

本文提出了一種計算非等徑顆粒視角系數(shù)的數(shù)值方法。該方法的關(guān)鍵是三種技術(shù)的新穎組合，即兩個非等徑顆粒視角系數(shù)的簡化積分表達(dá)式、球體上斐波那契積分方案和非均勻變換的組合。簡化積分不僅提供了一個簡單的積分公式，用于計算兩個非等徑顆粒之間的無阻礙視角系數(shù)的精確值，而且在提高計算有阻礙顆粒視角系數(shù)精度中起著基礎(chǔ)作用。斐波那契積分方案的張量積形式本質(zhì)上是一種光線跟蹤方法，對于計算視角系數(shù)的對偶積分而言，是一種計算量較少但具有足夠精度的方法。使用特定的非均勻變換可以很好地調(diào)節(jié)兩個積分的性質(zhì)，從而提高斐波那契積分法的求解精度。

對于具有阻礙顆粒的兩個非等徑顆粒的視角系數(shù)，從已經(jīng)實現(xiàn)的幾個算例評估了該方法的準(zhǔn)確性。此外，還以p1=1作為特例，將本文獲得的視角系數(shù)與文獻(xiàn)的結(jié)果進(jìn)行比較，顯示出了較好的一致性，尤其在NF=55的情況下。對于含有多顆粒粉末床中的視角系數(shù)，多個算例的結(jié)果表明，在有效區(qū)域尺寸D=6時，計算多顆非等徑顆粒視角系數(shù)的相對誤差為5%，并且與傳統(tǒng)的蒙特卡洛方法相比，計算效率提高約30%。