任意折射率介質(zhì)中的光束追跡方法

任意折射率介質(zhì)中的光束追跡方法

1網(wǎng)格節(jié)點(diǎn)不同時(shí)點(diǎn)折射波場(chǎng)模擬由于存在復(fù)雜的結(jié)構(gòu)，如激波、邊緣層、剪切層和水流，高速流動(dòng)場(chǎng)的流場(chǎng)介質(zhì)密度不均，沒(méi)有規(guī)律性，因此折射不可避免。計(jì)算流體力學(xué)所計(jì)算網(wǎng)格點(diǎn)或者網(wǎng)格單元中心點(diǎn)的密度(折射率)代表了該點(diǎn)處介質(zhì)的密度(折射率),所有網(wǎng)格點(diǎn)信息的集合則描述了整個(gè)流場(chǎng)的折射率分布,這是對(duì)無(wú)規(guī)則非均勻折射率場(chǎng)的一種描述方法。目前對(duì)折射率分布規(guī)律已知的介質(zhì)進(jìn)行光線追跡的研究已經(jīng)較為普遍,并且得到了廣泛的應(yīng)用,但是針對(duì)任意無(wú)規(guī)則折射率分布介質(zhì)的研究則相對(duì)較少。這些研究的基本思想是將每個(gè)網(wǎng)格單元視為等折射率單元,以兩個(gè)格子之間的界面作為折射面,根據(jù)Snell折射定律進(jìn)行光線追跡。然而在實(shí)際情況中,計(jì)算高速流場(chǎng)所采用的網(wǎng)格在大部分情況下是非均勻的,網(wǎng)格尺寸相差較大,而且存在激波等折射率變化劇烈的區(qū)域,不能簡(jiǎn)單認(rèn)為網(wǎng)格單元內(nèi)折射率處處相等,而且以一個(gè)網(wǎng)格單元作為一次追跡步長(zhǎng)也不合理。本文基于三維空間離散網(wǎng)格,網(wǎng)格節(jié)點(diǎn)介質(zhì)折射率已知,通過(guò)實(shí)時(shí)計(jì)算光傳播路徑上所需各點(diǎn)的折射率和折射率梯度,采用Runge-Kutta法實(shí)現(xiàn)光線追跡。在追跡過(guò)程中,每一步追跡的起點(diǎn)和終點(diǎn)不一定在網(wǎng)格點(diǎn)或者網(wǎng)格單元分界面上,并且追跡步長(zhǎng)根據(jù)當(dāng)?shù)卣凵渎侍荻群途W(wǎng)格尺寸自適應(yīng)調(diào)節(jié),以提高追跡的精度和效率。2限制節(jié)點(diǎn)法選取下式描述了折射率任意分布介質(zhì)中的光線傳播,是研究梯度折射率介質(zhì)光線傳播的基本方程之一:式中s為光線傳播路徑上的弧長(zhǎng),r為光線矢徑,n為折射率,▽n為折射率梯度。該方程在大多數(shù)情況下是難以求解的,僅對(duì)等折射率面為平面、球面和圓柱面等比較特殊的形式有解析解。對(duì)于任意非均勻折射率介質(zhì)中的光線傳播路徑,基本上沒(méi)有解析解,一般采用數(shù)值方法進(jìn)行求解。由于Runge-Kutta法在達(dá)到相同精度的情況下,花費(fèi)的運(yùn)算量較少,因此本文選用RungeKutta法作為基本光線追跡手段。對(duì)(1)式,令光線矢量Tray=ndr/ds=dr/dt,式中dt=ds/n是一個(gè)引入的參量,則光線傳播方程可寫為式中有3個(gè)分量,可用三元一維數(shù)列解得,記R為光線傳播路徑上的全坐標(biāo),T為光線矢量分量,D為(2)式右端項(xiàng)。經(jīng)過(guò)推導(dǎo)可得Runge-Kutta法計(jì)算公式如下:這里矩陣式中Δt=Δs/n,Δs為空間追跡步長(zhǎng)。當(dāng)已知光線初始條件R=R0(x0,y0,z0),T=T0時(shí),可以逐次進(jìn)行追跡,在每一步追跡完成后,需要計(jì)算該步末端點(diǎn)的折射率和折射率梯度,作為下一次追跡的起始條件,依此往后,直到完成整個(gè)光線追跡過(guò)程。3離散空間折射率及其梯度求解采用Runge-Kutta法對(duì)折射率空間任意離散分布的流場(chǎng)介質(zhì)進(jìn)行光線追跡時(shí),其追跡精度除了Runge-Kutta法本身精度的限制外,還受到3個(gè)因素的影響:1)追跡路徑上的折射率的計(jì)算精度;2)追跡路徑上的折射率梯度的計(jì)算精度;3)追跡步長(zhǎng)的大小。下面將分別對(duì)這3個(gè)因素展開(kāi)討論。3.1距離加權(quán)平均插值方法在離散空間內(nèi),已知網(wǎng)格節(jié)點(diǎn)上的折射率,需要求解其他空間任意點(diǎn)的折射率,通常采用插值方法。常用的插值方法有距離加權(quán)平均插值方法、補(bǔ)丁法、雙三次樣條插值和三線性插值法,此處選用精度和穩(wěn)定性均較好的距離加權(quán)平均插值方法。在三維離散空間內(nèi),如圖1所示網(wǎng)格單元,已知單元的八個(gè)頂點(diǎn)的折射率為ni(i=1,2,…,8),要求解單元內(nèi)部任一點(diǎn)P(x,y,z)的折射率nP。設(shè)點(diǎn)P到八個(gè)點(diǎn)(xi,yi,zi)的距離為di,距離加權(quán)平均插值方法可表示為3.2充分利用有限差分方法求解流場(chǎng)已知空間網(wǎng)格節(jié)點(diǎn)的折射率,求解網(wǎng)格內(nèi)任意點(diǎn)的折射率梯度,需要分兩步:第一步求解網(wǎng)格點(diǎn)本身的梯度值,第二步利用網(wǎng)格單元頂點(diǎn)的梯度插值求解網(wǎng)格內(nèi)任意點(diǎn)的梯度值。其中第二步的插值方法與求解折射率的方法相同,這里重點(diǎn)論述第一步的求解方法。求解網(wǎng)格節(jié)點(diǎn)的梯度需要利用節(jié)點(diǎn)周圍其它多個(gè)點(diǎn)進(jìn)行求解,常用的方法有Horn梯度算子,Barron梯度算子等,此處介紹計(jì)算流體力學(xué)中具有二階精度的有限差分方法求解梯度。計(jì)算流體力學(xué)一般采用曲線坐標(biāo)系(ξ,η,ζ)下的N-S方程求解流場(chǎng),如圖2所示,而光線追跡一般在笛卡兒坐標(biāo)系(x,y,z)下進(jìn)行,所以需要求解笛卡兒坐標(biāo)下的折射率梯度▽n。在曲線坐標(biāo)系下,首先運(yùn)用有限差分方法求得折射率梯度▽n的三個(gè)分量▽nξ,▽nη,▽nζ,再利用Jacobi坐標(biāo)變換矩陣得到笛卡兒坐標(biāo)下的分量▽nx,▽ny,▽nz。1)求解▽nξ,▽nη,▽nζ以ξ方向?yàn)槔?采用如下的差分形式求解網(wǎng)格點(diǎn)(i,j,k)處的折射率的梯度左邊界網(wǎng)格點(diǎn)采用三點(diǎn)向前偏心差分中間網(wǎng)格點(diǎn)采用中心差分右邊界網(wǎng)格點(diǎn)采用三點(diǎn)向后偏心差分上面各式中,下標(biāo)i,j,k與ξ,η,ζ的關(guān)系如圖2所示。2)利用Jacobi坐標(biāo)變換矩陣得到▽nx,▽ny,▽nz在三維空間中,Jacobi坐標(biāo)變換矩陣為利用Jacobi矩陣即可求得折射率梯度在笛卡兒坐標(biāo)x,y,z三個(gè)方向上的分量:3.3追跡步長(zhǎng)逐步步長(zhǎng)為小地區(qū),追跡步長(zhǎng)根據(jù)實(shí)際追跡算法合理在光線追跡過(guò)程中,一般給定追跡步長(zhǎng)Δs,Δs越小則追跡的精度越高,但這將增加計(jì)算量,所以在實(shí)際追跡過(guò)程中可考慮對(duì)追跡步長(zhǎng)進(jìn)行自適應(yīng)調(diào)整,在折射率變化大的地方,步長(zhǎng)小;在折射率變化小的地方,步長(zhǎng)可以大一些,這樣既保證精度又提高了效率。3.3.1步長(zhǎng)調(diào)節(jié)函數(shù)與雙梯度梯度轉(zhuǎn)化改變追跡空間步長(zhǎng)的一種方法是給定不變的光學(xué)步長(zhǎng)Δv=nΔs,根據(jù)當(dāng)?shù)卣凵渎蕁(x,y,z)變化調(diào)節(jié)追跡空間步長(zhǎng)Δs。但這種方法存在不足:當(dāng)折射率n本身較大而又基本不變化,就會(huì)導(dǎo)致空間步長(zhǎng)Δs較小,對(duì)精度提高并沒(méi)有明顯改善,反而增加了計(jì)算時(shí)間。本質(zhì)上,光線對(duì)折射率梯度▽n最為敏感,為了提高光線追跡的精度,應(yīng)該重點(diǎn)考慮折射率梯度的影響,文獻(xiàn)的方法也體現(xiàn)了這一思想。本文基于相同的出發(fā)點(diǎn),即在折射率變化比較劇烈(梯度大)的地方減小追跡步長(zhǎng)對(duì)光線進(jìn)行細(xì)致追跡,但采用不同的方式實(shí)現(xiàn),無(wú)需對(duì)網(wǎng)格單元進(jìn)行細(xì)分和重構(gòu)。其方法是利用網(wǎng)格折射率梯度xg作為變量,構(gòu)造步長(zhǎng)調(diào)節(jié)函數(shù)f(xg)。由于調(diào)節(jié)函數(shù)難以給出解析表達(dá)式,經(jīng)過(guò)大量的數(shù)值計(jì)算實(shí)驗(yàn),給出如下分段函數(shù)形式的步長(zhǎng)調(diào)節(jié)函數(shù)追跡步長(zhǎng)Δs=f(xg),Δi=(ai×bi×ci)/3為當(dāng)?shù)鼐W(wǎng)格單元i的平均幾何尺寸,ai,bi,ci為當(dāng)?shù)鼐W(wǎng)格單元邊長(zhǎng),如前面圖2所示。對(duì)(11)式有三點(diǎn)說(shuō)明:1)網(wǎng)格折射率梯度xg是指相近網(wǎng)格單元之間的折射率變化幅度(單位:單元格),不是當(dāng)?shù)氐慕^對(duì)梯度值(單位:m);2)Δi為追跡步長(zhǎng)的最大值,即使在折射率梯度很小的情況下,限定每次追跡步長(zhǎng)不超過(guò)當(dāng)?shù)鼐W(wǎng)格平均幾何尺寸;3)最小的追跡步長(zhǎng)是當(dāng)?shù)鼐W(wǎng)格平均幾何尺寸Δi的0.05倍,這是對(duì)存在激波這種極端情況進(jìn)行數(shù)值考察后給出的參考數(shù)值,限于篇幅,此處不再給出詳細(xì)考察過(guò)程。3.3.2單元內(nèi)分辨率q由于采用了變步長(zhǎng)追跡,每一步追跡完成后的位置并不一定在網(wǎng)格點(diǎn)或者網(wǎng)格單元分界面上,需要界定每一步光線末端點(diǎn)所在的網(wǎng)格單元,據(jù)此插值求解該點(diǎn)的折射率和折射率梯度,方法如下:判斷點(diǎn)P是否在某個(gè)網(wǎng)格單元K內(nèi),可通過(guò)計(jì)算該點(diǎn)與單元K的八個(gè)頂點(diǎn)順次組成的12個(gè)四面體體積之和VP與該單元體的體積VK的關(guān)系進(jìn)行判別。考慮實(shí)際計(jì)算存在誤差,需要給定一個(gè)誤差限ε進(jìn)行判斷:4x-y-z方向的光學(xué)方向余弦在下面計(jì)算中,以z軸為光軸,在光軸坐標(biāo)zi處對(duì)應(yīng)的光線真實(shí)傳播軌跡的矢徑為ri,經(jīng)過(guò)若干步追跡后的計(jì)算值為,定義計(jì)算值與精確解的相對(duì)誤差考察標(biāo)準(zhǔn),折射率徑向變化規(guī)律如下:在(14)式對(duì)應(yīng)折射率分布的介質(zhì)中,光線傳播軌跡的精確解為式中x0,y0,z0為初始點(diǎn)的坐標(biāo),p0,q0,L0為初始點(diǎn)對(duì)應(yīng)x,y,z方向的光學(xué)方向余弦,算例中所涉及到的坐標(biāo)長(zhǎng)度單位均以mm度量。在給定折射率空間離散分布時(shí),考慮到折射率為徑向分布,x,y方向網(wǎng)格間距Δx,Δy需要小一些,z方向網(wǎng)格間距Δz可以大一些。給定三個(gè)空間離散網(wǎng)格形式:1)網(wǎng)格1:△x=△y=1mm,Δz=5mm;2)網(wǎng)格2:Δx=Δy=2mm,Δz=10mm;3)網(wǎng)格3:Δx=Δy=4mm,Δz=50mm。為了充分考察本光線追跡方法的精度,此處沿z方向10m范圍進(jìn)行追跡計(jì)算。不失一般性,假設(shè)光線初始入射點(diǎn)坐標(biāo)為(1.0,1.0,0),取n(0)=1.0,光線初始點(diǎn)光學(xué)方向余弦為(0,0,1.0),考慮不同追跡步長(zhǎng),計(jì)算得到的結(jié)果如圖3所示。4.1網(wǎng)格追跡精度選擇折射率徑向分布時(shí)光線傳播軌跡精確解為由圖3可知,在同一個(gè)網(wǎng)格下,對(duì)于追跡步長(zhǎng)分別為0.2mm,1.0mm以及步長(zhǎng)自適應(yīng)追跡三種情況,其追跡精度保持一致,差別不明顯;隨著追跡步長(zhǎng)增大所需的步數(shù)減少,步長(zhǎng)自適應(yīng)方法的追跡步數(shù)最少,效率最高。此現(xiàn)象表明,光線追跡的精度在已確定折射率空間分布情況(即同一網(wǎng)格情況)下,對(duì)追跡步長(zhǎng)不是很敏感,但追跡的效率隨步長(zhǎng)的增大而提高。對(duì)于不同的網(wǎng)格,追跡精度變化則比較明顯,三個(gè)網(wǎng)格追跡的精度為10-5～10-4量級(jí),如圖4所示。這主要是因?yàn)?稀疏的網(wǎng)格節(jié)點(diǎn)信息對(duì)折射率場(chǎng)的描述不精細(xì),在插值求解折射率和折射率梯度時(shí)造成的誤差較大,而追跡方法本身造成的誤差是次要原因。在達(dá)到同樣精度的條件下,隨著網(wǎng)格變得稀疏,追跡的步數(shù)明顯降低,效率顯著提高;而且追跡步長(zhǎng)由空間網(wǎng)格尺寸和網(wǎng)格折射率梯度得到,無(wú)需事先評(píng)估追跡步長(zhǎng)到底給定多少才能滿足精度要求。另外,折射率徑向分布中的參數(shù)α反映了折射率的變化幅度(即梯度)?；诰W(wǎng)格1(gird1),針對(duì)不同的α值,考察本方法對(duì)不同折射率變化幅度時(shí)的追跡精度如圖5所示。為了直觀表示,圖5縱坐標(biāo)采用對(duì)數(shù)坐標(biāo)。結(jié)果表明,在同一個(gè)網(wǎng)格下(即折射率場(chǎng)空間分布一定),折射率變化幅度參數(shù)α每增加一個(gè)量級(jí)(從10-6～10-4),折射率梯度依次增加二個(gè)量級(jí)(從10-12～10-8);沿z方向追跡10m后,折射率變化Δn的范圍分別達(dá)到10-7,10-10與10-12量級(jí),追跡精度依次降低二個(gè)量級(jí)(從10-9～10-5)。根據(jù)前面分析可知,由于折射率空間離散網(wǎng)格沒(méi)有隨折射率變化劇烈而加密,導(dǎo)致折射率場(chǎng)信息不精細(xì),從而導(dǎo)致追跡精度降低。4.2與模擬傳播路徑比較選取螺旋光線為對(duì)象進(jìn)行追蹤,對(duì)應(yīng)的折射率分布為對(duì)應(yīng)的螺旋光線軌跡的精確解為初始條件p0=0,。同樣基于網(wǎng)格1,取α=0.01,n(0)=1.5,將計(jì)算得到的結(jié)果與光線真實(shí)傳播路徑進(jìn)行比較,如圖6所示,兩者吻合的程度比較高。圖7為追跡過(guò)程中對(duì)應(yīng)誤差,在沿z軸追跡10m計(jì)算過(guò)程中,誤差始終保持在10-4量級(jí)。應(yīng)該說(shuō)對(duì)于螺旋光線這類特殊光線的追跡效果是較好的,由此說(shuō)明本方法具有較高的精度。5有限差分追跡本文所建立的光線追跡方法是針對(duì)折射率任意空間離散分布介質(zhì)中的光線追跡,可以直接應(yīng)用于計(jì)算流體力學(xué)的網(wǎng)格,無(wú)需再重新劃分和重構(gòu)追跡所需的網(wǎng)格;網(wǎng)格節(jié)點(diǎn)上的折射率梯度采用計(jì)算流體力學(xué)中具有二階精度的有限差分方法進(jìn)行求解,追跡過(guò)程中所需折射率和折射率梯度采用距離加權(quán)插值得到,同時(shí)提