液體激光器熱流場的基本理論

一、熱傳遞及流體力學原理：熱量傳遞由三種方式：熱傳導、熱對流和熱輻射。1.熱傳導熱量有物體的一部分傳到另一部分，或由一個物體傳到與其接觸的另一個物體的傳熱現(xiàn)象叫熱傳導。其特點是，熱量流動而物體各部分仍保持著宏觀的靜止。熱傳導定律——傅里葉假設傅里葉假設：在單位時間內通過微元等溫面dA的熱量dQ與成正比，即dQ二_dQ二_k—ndA

dn(1-1)式中：“-”號——熱量是從溫度高的地方流向溫度低的地方，其方向與溫度梯度相反；傳熱系數(shù)，它是表征材料導熱能力的物理參數(shù)，單位是W/(m2?K)；冬 溫度場梯度值；dnn——等溫面法線方向的單位矢量。熱流強度單位時間通過單位等溫面積的熱量叫熱流強度。按照傅里葉假設，熱傳導熱流強度的計算公式為dTq=dQ/dA=一k——n (1-2)dn式中q為熱傳導熱流強度，單位是W/m2。熱傳導熱流強度q是一個矢量。2.對流換熱通過流體的流動而換熱的現(xiàn)象叫對流換熱。對流換熱是一個復雜的過程，它既包括流體的對流作用，也包含流體的導熱作用。按牛頓公式對流換熱的熱流強度的計算公式為q=h(T_T)n (1-3)ces式中：q——對流熱流強度，單位是W/m2；ch——表面?zhèn)鳠嵯禂?shù)，單位是W/(m2?K)；T——流體的平均溫度；eT——固體表面的溫度；sn——指向固體表面的單位法向矢量。對流熱流強度q是一個矢量，其方向垂直于固體表面，以指向固體為c正。對流系數(shù)h是考慮了影響對熱換熱的各種因素的一個綜合系數(shù)，這些因素是流體的速度、流體的物理特征、換熱表面的形狀和尺寸等。熱輻射在本研究中不做考慮。3.管內流體對流換熱a.對流速度邊界層當具有粘性的流體流過壁面時，就會在壁面上產生粘滯力。粘滯力阻礙了流體的運動，使靠近壁面的速度降低，使直接貼附于壁面的流體近于停滯不動。如圖1所示，沿壁面發(fā)現(xiàn)y方向不同點的速度v，就是速度分布圖。速度分布圖表明：從y=0處v=0開始，速度v隨著y方向的增加而迅速增大;當經(jīng)過厚度為8的薄層時，v就接近主流速度v。把所有y=8的點連接成線即為速度邊界層，8為邊界層厚度。邊界層圖1.速度邊界層b.層流、湍流流體進入管口一段距離后，管壁兩側的邊界層才能在管中心回合，這時管斷面流速分布和流動狀態(tài)才達到定型，這段距離通常稱為入口段。之后的流速分布不再改變，流態(tài)定型，流動達到充分發(fā)展，稱為流動充分段。如圖2、圖3所示。

圖2層流狀態(tài)的管內流動流速分布v入口段圖2層流狀態(tài)的管內流動流速分布v入口段充分發(fā)展段(1-4)圖3紊流狀態(tài)的管內流動流速分布層流和紊流的流態(tài)判斷由雷諾準則（1-4）式判定ud(1-4)Re=v式中 Re——雷諾數(shù)；v 運動粘度，單位m2Js；u——流體速度；d 圓管內徑（方管則約為截面面積丘x2）,單位m當Re<2320時為層流；當2320vRevl04時為過渡時期；當104<Re時為旺盛紊流。層流時，速度分布曲線為拋物線形狀；紊流時，速度分布曲線呈對數(shù)曲線形狀。管子入口段的長度、彎曲管道對流場的影響入口段管內流動的流速分布不穩(wěn)定，并且隨著流距而改變；當流體在彎管內流動時，因離心力而產生二次環(huán)流，增加了流動中的擾動，使流動場變?yōu)榉浅浞职l(fā)展流動，并且彎曲半徑越小，影響越大。介于上面的情況，如果在分析中希望達到充分發(fā)展以保持流場的穩(wěn)定，必須增長管的長度（彎曲后管的長度）。在紊流狀態(tài)下管長與管內徑之比->50時（層流要小些），可以認為d管內流體已經(jīng)達到充分發(fā)展，即為穩(wěn)定流場。二、有限元1.原理將一個表示結構或連續(xù)體的求解域離散為若干個單元，并通過它們的邊界上的節(jié)點相互聯(lián)結成為組合體。用每個單元內所假設的近似函數(shù)來分片地表示全求解域內待求的未知場變量。而每個單元內的近似函數(shù)由未知場函數(shù)(或及其導數(shù)，為敘述方便，后面略去此加注)在單元各個結點上的數(shù)值和與其對應的插值函數(shù)來表達。由于聯(lián)結相鄰單元的結點上，場函數(shù)應具有相同的數(shù)值，銀耳將它們用作數(shù)值求解的基本未知量。這樣一來，求解原來待求場函數(shù)的無窮多自由度問題轉換為求解場函數(shù)結點值的有限自由度問題。通過和原問題數(shù)學模型(基本方程、邊界條件)等效的辯分原理或加權余量法，建立求解基本未知量(場函數(shù)的結點值)的代數(shù)方程組或常微分方程組。2.有限元解決熱傳導問題a.引言在一般的三維問題中，瞬態(tài)溫度場的場變量0(x,y,z,)在直角坐標中應滿足的微分方程是Pcdt dxPcdt dxvxdx丿dydzvzQz-PQ二0(在o內)(2-1)邊界條件是Q0 邊界條件是Q0 Q0 Q0kn+kn+kn=qxQxx yQyy zQzz(在r邊界上) (2-2)1(在r邊界上) (2-3)2Q0 Q0 Q0 Q0 Q0kn+kn+k nxQxx yQyyzQzz二h(0-0)a(在r邊界上)3(2-4)式中：P 材料的密度，單位是kg/m3；c 材料的比熱容，單位是J/(kg-K)；0=0(r,t)——在r邊界上的給定溫度；1q=q(r,t)——在r邊界上的給定熱流密度，單位是w/m2；20=0(r,t)——對于r邊界，在自然對流條件下，0是外界環(huán)境溫度;aa3a在強迫對流條件下，0是邊界層的絕熱壁溫度。a微分方程(2-1)是熱量平衡方程。式中第1項是微體升溫需要的熱量；第2,3,4項是由x,y和z方向傳入微體的熱量；最后一項是微體內熱源產生的

熱量。微分方程表明：微體升溫所需的熱量應與傳入微體的熱量以及微體內熱源產生的熱量相平衡。(2-2)式是在r邊界上給定溫度O=e(r,t),稱為第一類邊界條件，它是1強制邊界條件。(2-3)式是在r邊界上給定熱流量q=q(r,t),稱為第二類邊2界條件，當q=0時就是絕熱邊界條件。(2-4)式是在r邊界上給定對流換熱3的條件，稱為第三類邊界條件。第二、三類邊界條件是自然邊界條件。r+r+r二r，r是域o的全部邊界。123當在一個方向上，例如z方向溫度變化為零時，方程(2-1)就退化為二維問題的熱傳導微分方程：pc也

apc也

ata匚a?)——k——axvxax丿k汨-PQ=0Vyay丿(在o內) (2-5)這時當場變量?(x,y,t)不再是z的函數(shù)，場變量應滿足的邊界條件是:帖?(r,t)+k竺n帖?(r,t)+k竺n=q(r,t)yayyk越nxax+k越n=h(?-?)yayya(在r邊界上)1(在r邊界上)2(在r邊界上)3在(2-5)~(2-8)式中，各項符號意義與在(2-1)~(2-4)式中的相同。求解瞬態(tài)溫度場問題是求解在初始條件下，即在?(x,y,z,0)=?0(x,y,z(2-6)(2-7)(2-8)(2-9)條件下滿足瞬態(tài)熱傳導方程及邊界條件的場函數(shù)?，?應是坐標和時間的函數(shù)。b.瞬態(tài)熱傳導問題1)瞬態(tài)熱傳導有限元的一般格式首先建立三維瞬態(tài)熱傳導問題的微分方程(2-1)式和邊界條件(2-2)，(2-3)和(2-4)式的等效積分形式，即

3e3c3e] 3匚3ewpc——-——Ik——-——k——3e3c3e] 3匚3ewpc——-——Ik——-——k——9t 9x(x3x丿9y帥JwG-e)dT+JT11JwT3 3k理nx3xx3(帥)_k丄-pQdQ+3zIz3z丿 _7r73e73ek——n+k——n+k——n-qdr+x3xxy3yyz3zz丿wr22+k遊n+k遊n-h(e-0)]dr=0y3yyz3zza丿(2-10)其中w，w，12 3 1已滿足?=&，則w=0，并且不失一般地可令1(2-11)W=W=W=8e(2-11)238efPc越〕+98e帥)98^(帥)38^k—+ k8efPc越〕+98e帥)98^(帥)38^k—+ k—+ k3e]-80pQdQ3z(z3z丿JQ+J8eqdr-J8eh(e-e)dr=o廠 廠 a丄2 丄3利用(2-12)式可以建立瞬態(tài)溫度場有限元的一般格式。首先將空間域Q9t丿9x(x9x丿k——+3yIy9y丿(2-12)離散為有限個單元體，在典型單元內溫度e仍可以近似地用結點溫度e插值i得到，但是注意此時的結點溫度是時間的函數(shù)，即e=e=*n(x,y,z)e(t) (2-13)iii=1插值函數(shù)N只是空間域的函數(shù)，它應具有插值函數(shù)的基本性質。將上式i代入(2-12)式，并且考慮8e的任意性，就可以得到用來確定n個結點溫度eii的有限元求解方程ce+ke=p (2-14)這是一組以時間t以獨立變量的線性常微分方程組。式中：C——熱容矩陣；K——熱傳導矩陣，C和K(在引入給定溫度以后)都是對稱正定矩陣；P——溫度在和列陣；e——結點溫度列陣；ef=業(yè)］——結點溫度對時間的導數(shù)列陣。Idt丿

矩陣C,K和P的元素由單元的相應的矩陣元素集成，即K丄Ke+ZHeijeC=工ije矩陣C,K和P的元素由單元的相應的矩陣元素集成，即K丄Ke+ZHeijeC=工ijeKe+He

ij ijeCeij(2-15)式中：P=工iPe+QieKeijPe+ Peqi Hiee單元對熱傳導矩陣的貢獻Heij單元熱交換邊界對熱傳導矩陣的修改；Ceij單元對熱容矩陣的貢獻；PeQi單元熱源產生的溫度載荷；Pe

qi單元給定熱流邊界的溫度載荷；廠6N6N 6N6N 6N6N)Ke=fij3(2-16)He=ihNNdrijreij3Ce=fpcNNdoij oe ijPe=fpQNdoQioe iPe=fqNdrqKe=fij3(2-16)He=ihNNdrijreij3Ce=fpcNNdoij oe ijPe=fpQNdoQioe iPe=fqNdrqi re2iPe=fhNdrH reaii re3(2-17)(2-18)(2-19)(2-20)(2-21)至此，已經(jīng)將時間域和空間域的偏微分方程問題在空間域內離散為N個結點溫度0(t)的常微分方程的初值問題。對于給定溫度值的邊界r上n個結i 1 1點，方程組(2-14)式中相應項應引入的條件是(2-22)0=0(i=1,2,…,n)(2-22)ii 1式中，i是r上n個結點的編號。112)一階常微分方程的求解瞬態(tài)熱傳導問題的有限元求解方程(2-14)式是一階常微分方程組。以下討論有限元的直接積分法求解該方程。直接積分法是指求解常微分方程組前不對它的形式進行變換，而是直接對它進行數(shù)值積分。它通?；趦蓚€基本概念，一是求解的時間域0<t<T劃分為若干個時間步長At，在一定數(shù)目的At時間區(qū)域內，假設Q和$的函數(shù)形式來近似方程的精確解；二是以僅在相隔At的離散時間點上滿足微分方程來代替時間域內任何時刻t都滿足微分方程。在以下的討論中首先將時間域0?T等分(不限于等分，這里是為了討論方便)為M個時間步長，At二TIM，并認為t=0時的初始溫度列陣$是已知0的。進一步假設t=0，t(二At)，t(二2At),???,t(=ntA)時刻的解$已經(jīng)TOC\o"1-5"\h\z1 2 n n求得。下一步要計算的是t(=t+At)時刻的$ 。從該求解過程建立起求解n+1 n n+1所有離散時間點場函數(shù)的一般算法步驟。2.1)用加權余量法建立兩點循環(huán)公式對于只有一階導數(shù)的常微分方程組(2-14)式，可以再兩個時間點t和tn n+1之間的At時間區(qū)域內，假設$和$采取如圖2所示的線性插值形式，即$(t+S)At)=N$+N$ (2-23)n nn n+1n+1其中N=1—爲 (0<g<1)nN=￡n+1(2-23)式對時間t求導，可得(2-24)$=N$+N$(2-24)nn n+1n+1其中1Atn+n+1AtJ(a)(b)(c)o=2(d)(e)圖1建立兩點循環(huán)公式的插值函數(shù)及權函數(shù)(a)(b)(c)o=2(d)(e)圖1建立兩點循環(huán)公式的插值函數(shù)及權函數(shù)由于采用(2-23)式的近似插值，在時間區(qū)域At內，方程(2-14)式將產生余量。對于這一時間區(qū)域，典型的加權余量格式可以表示為如下形式：n ，"‘丄 n ，"‘丄 n+1丿+kN?+N?n n+1)n ，"‘丄 n+1丿-Pdg二0(2-25)當求解初值問題時，如果已知一組參數(shù)?，則可以利用(2-25)式近似確n定另一組參數(shù)?。將(2-23)和(2-24)式代入(2-25)式就可以得到時間區(qū)域n+1At前后結點上兩組參量的關系式，即fKJ1fKJ1誕dg+CJ1w朮'Vo oAt丿-J1wPd^=000+fKJ1w(1弋)dg-CJ1w朮'00At丿n+1(2-26)式中可以代入不同的權函數(shù)。在以上討論中假定熱傳導矩陣K和熱容矩陣C不隨時間t而變化。(2-26)式可以表示為任何權函數(shù)都使用的一般形式，即(C/At+K0)0 +「-C/At+K(1-8)0]=P (2-27)n+1 n式中(2-28)J1誕d^/f1wdgJ1wPd^.J1wdg0■'0(2-28)一種方便的做法就是假定P采用與未知場函數(shù)0相同的插值表達式，這時將得到P=P0+P(1-0) (2-29)TOC\o"1-5"\h\zn+1 n當0和P都已知時，就可以由(2-27)式求得下一時刻的0 。這就是連點循環(huán)nn n+1公式。可以更明顯地表示成(2-30)(2-31)K0=Q (n=0,1,2,M)(2-30)(2-31)\o"CurrentDocument"n+1 n+1其中K=C/At+0KQ=「C/At-(l-0)K]0+(1-0)P+0P>n+1「 」n n n+1丿利用上式，從t=0出發(fā)，提一次一次遞推求得結點溫度列陣0(t)的各個瞬時值0,0，…,0。1 2M2.2)算法步驟利用直接積分的兩點循環(huán)公式遞推求解瞬態(tài)熱傳導問題的有限元微分方程的算法步驟可以歸結如下：初始計算Q形成系統(tǒng)系數(shù)矩陣C和K；Q給定0；0Q選擇參數(shù)0和時間步長At；Q形成系統(tǒng)的有效系數(shù)矩陣K=C/At+0K；

⑨三角分解K,K=LDLT.對于每一時間步長：Q形成向量P；n+19形成有效向量Qn+1Pn+1二「C/At-(1-O)k伸+(i-e)9形成有效向量Qn+1Pn+1TOC\o"1-5"\h\zn n\o"CurrentDocument"9回代求解0 ,LDL0 =Q。n+1 n+1 n+12.3)參數(shù)0的選擇以上通過加權余量法得到的求解瞬態(tài)熱傳導方程的兩點循環(huán)公式(2-27)本質上是一組加權的差分公式，因為在每個時間區(qū)域At內，對于(t+0At)點n(0<0<1)建立的差分公式為(2-32)0(t+0At)=(l-0)0+00(2-32)n n n+10(t+0At)=(0n n+1(2-33)該兩式實際上就是利用時間區(qū)域At內線性插值公式得到的(2-23)和(2-24)式,只是坐標參數(shù)g換成0而已。將(2-32)和(2-33)式代入(2-14),并將P表示成和0相同的差分形式，則得到建立于(t+0At)時刻的差分方程。它和前面利用加權余量法的兩點循環(huán)公n式(2-27)式相同。從此可以清楚地認識0的物理意義，0的取值決定了在At時間區(qū)域內建立差分方程的具體地點。從圖2給出的一組權函數(shù)及相應的0值1的關系也可以看出0的物理意義。圖2中的(a)、(b)、(c)集中在點n，n+丄，2n+1上加權，得到的是著名的前差分公式(歐拉差分公式)，中心差分公式(Crank-Nicholson差分公式)和后差分公式。(d)是在時間區(qū)域內權函數(shù)等于常數(shù)，其結果和中心差分的結果相同。(e)和⑴為伽遼金型的權函數(shù)，其結果分12別和權函