本發(fā)明涉及半導(dǎo)體溫差發(fā)電技術(shù)領(lǐng)域，尤其涉及一種熱端溫度時(shí)變的動(dòng)態(tài)溫差發(fā)電系統(tǒng)時(shí)域分析方法。

背景技術(shù)：

1821年德國(guó)科學(xué)家seebeck發(fā)現(xiàn)，在兩種不同的金屬構(gòu)成的閉合回路中，當(dāng)兩個(gè)接頭存在溫差時(shí)，回路將產(chǎn)生電流，這個(gè)現(xiàn)象被稱為塞貝克效應(yīng)。溫差發(fā)電就是利用塞貝克效應(yīng)，通過在熱電材料兩端維持一定的溫差，從而產(chǎn)生一定的電壓和電功率輸出。通過研究發(fā)現(xiàn)，半導(dǎo)體材料的熱電優(yōu)值較大，目前人們?cè)跍夭钕到y(tǒng)中應(yīng)用的熱電材料都是半導(dǎo)體材料，所以又稱為半導(dǎo)體溫差發(fā)電。

隨著人們對(duì)能源危機(jī)的認(rèn)識(shí)，溫差發(fā)電技術(shù)可以利用自然界中存在的大量的溫差以及工業(yè)余熱，具有良好的綜合社會(huì)經(jīng)濟(jì)效益。同時(shí)，隨著人們對(duì)空間探索興趣，醫(yī)用物理學(xué)的進(jìn)展，和大規(guī)模無線傳感器的應(yīng)用，需要開發(fā)一類能夠自身供給能量并無需照看的電源系統(tǒng)，顯然，溫差發(fā)電對(duì)這些應(yīng)用極為合適。溫差發(fā)電作為一種全固態(tài)能量轉(zhuǎn)換方式，具有無介質(zhì)泄露，無噪聲，性能可靠，維護(hù)少等優(yōu)點(diǎn)的綠色環(huán)保能源，在微型能源、低品位能源、廢能源利用方面的應(yīng)用價(jià)值越來越明顯。因而，盡快實(shí)現(xiàn)溫差發(fā)電技術(shù)及其應(yīng)用的產(chǎn)業(yè)化具有重要的現(xiàn)實(shí)意義。而溫差發(fā)電系統(tǒng)的分析和設(shè)計(jì)依賴于實(shí)際的工作狀況，故針對(duì)不同的實(shí)際情況進(jìn)行合理的溫差發(fā)電系統(tǒng)模型是對(duì)溫差電系統(tǒng)的各種參數(shù)進(jìn)行更具體的分析和設(shè)計(jì)的基礎(chǔ)。

中國(guó)專利文獻(xiàn)cn201410189306.6公開了一種溫差發(fā)電系統(tǒng)電動(dòng)勢(shì)計(jì)算方法，在該專利中采用的對(duì)溫差電系統(tǒng)進(jìn)行分析計(jì)算的數(shù)學(xué)模型，對(duì)實(shí)際溫差發(fā)電系統(tǒng)進(jìn)行了如下的近似：

認(rèn)為溫差電系統(tǒng)中的物理特性參數(shù)是恒定值，即為常數(shù)；

認(rèn)為邊界處的熱流或溫度是時(shí)不變的，即為常數(shù)，

故專利文獻(xiàn)cn201410189306.6中溫差電系統(tǒng)電動(dòng)勢(shì)的計(jì)算方法，將不同物質(zhì)交界處的邊界條件，簡(jiǎn)單近似為連接處的流體平衡，只適用于溫差發(fā)電系統(tǒng)是處于常物性的，邊界處可簡(jiǎn)單近似為平衡態(tài)情況下，穩(wěn)態(tài)線性溫差電系統(tǒng)的分析，而對(duì)于溫差發(fā)電系統(tǒng)處于非平衡態(tài)時(shí)是不適用的，即不能用于分析動(dòng)態(tài)溫差發(fā)電系統(tǒng)的瞬態(tài)過程。

而實(shí)際的溫差電系統(tǒng)中，如果根據(jù)工況條件，溫差發(fā)電系統(tǒng)的邊界處熱源需等效為隨時(shí)間變化的溫度輸入時(shí)，此時(shí)的溫差發(fā)電系統(tǒng)處于非平衡態(tài)，溫差發(fā)電系統(tǒng)中不同物質(zhì)交界處的邊界條件不能簡(jiǎn)單近似為連接處的流體平衡。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是提供一種熱端溫度時(shí)變的動(dòng)態(tài)溫差發(fā)電系統(tǒng)時(shí)域分析方法，通過本發(fā)明中的方法得到動(dòng)態(tài)溫差發(fā)電內(nèi)部溫度和電場(chǎng)強(qiáng)度的數(shù)值解，可進(jìn)一步對(duì)熱端熱流時(shí)變動(dòng)態(tài)溫差發(fā)電系統(tǒng)的進(jìn)行更具體的分析和設(shè)計(jì)，并且結(jié)果更加精確。

本發(fā)明采用的技術(shù)方案：

一種熱端溫度時(shí)變的動(dòng)態(tài)溫差發(fā)電系統(tǒng)時(shí)域分析方法，具體包括以下幾個(gè)步驟：

步驟一，根據(jù)溫差發(fā)電系統(tǒng)材料的基本物理特性，熱力學(xué)和電學(xué)基本定律，導(dǎo)出由d對(duì)半導(dǎo)體溫差熱電偶構(gòu)成的溫差發(fā)電系統(tǒng)熱路和電路的數(shù)學(xué)模型；

熱路中各個(gè)區(qū)域的數(shù)學(xué)模型為：

區(qū)域i中，

區(qū)域ii中，

區(qū)域iii中，

電路各個(gè)區(qū)域的數(shù)學(xué)模型為：

區(qū)域ii中，

區(qū)域iv中，

其中，溫差電系統(tǒng)熱路分成三個(gè)區(qū)域：區(qū)域i，區(qū)域ii和區(qū)域iii，區(qū)域i代表熱端散熱器，區(qū)域ii代表d對(duì)半導(dǎo)體溫差熱電偶，區(qū)域iii代表冷端散熱器；溫差電系統(tǒng)中電路分成兩個(gè)區(qū)域：區(qū)域ii和區(qū)域iv；區(qū)域ii代表d對(duì)半導(dǎo)體溫差熱電偶，區(qū)域iv代表負(fù)載電阻rl；

其中，ρi和ρiii分別為區(qū)域i和iii物質(zhì)的密度，cvi和cviii分別為區(qū)域i和iii物質(zhì)的定容比熱，ki和kiii分別為區(qū)域i和iii物質(zhì)的熱導(dǎo)率。t為溫度，為溫度梯度；e為電場(chǎng)強(qiáng)度，ρii分別為區(qū)域ii物質(zhì)的密度，cvii為區(qū)域ii物質(zhì)的定容比熱，εii為區(qū)域ii物質(zhì)的介電常數(shù)，j為電路中的電流密度，電流j0為t＝0時(shí)刻電路中的的電流密度，αii，σii，kii，分別為區(qū)域ii物質(zhì)的總等效塞貝克系數(shù)，總等效電導(dǎo)率，總等效熱導(dǎo)率；αii＝d*(αp-αn)，kii＝m*(kp+kn)；αp，σp，kp分別為溫差熱電偶中p型半導(dǎo)體材料的塞貝克系數(shù)，電導(dǎo)率和熱導(dǎo)率，αn，σn，kn，分別為溫差熱電偶中n型半導(dǎo)體材料的塞貝克系數(shù)，電導(dǎo)率和熱導(dǎo)率；σiv分別為負(fù)載電阻的電導(dǎo)率。

步驟二，根據(jù)實(shí)際工況條件確定熱端溫度時(shí)變的動(dòng)態(tài)溫差發(fā)電系統(tǒng)的熱路和電路的邊界條件：

熱路中的邊界條件為：

邊界a:根據(jù)熱源提供能量的形式，熱端等效為溫度為th(t)；

邊界b:通過邊界處的熱流密度為連續(xù)的，溫度為連續(xù)的；

邊界c:通過邊界處的熱流密度為連續(xù)的，溫度為連續(xù)的；

邊界d：根據(jù)冷源端散熱的形式，冷端等效為溫度為tl(t)；

電路中的邊界條件為：

邊界m:通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的；

邊界n:通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的；

其中熱路中區(qū)域i和熱源的交界為邊界a，區(qū)域i和區(qū)域ii的交界為邊界b，區(qū)域ii和區(qū)域iii的交界為邊界c，區(qū)域iii和冷源的交界為邊界d；其中電路中區(qū)域ii在熱源端和區(qū)域iv的交界為邊界m，其中區(qū)域ii在冷源端和區(qū)域iv的交界為邊界n；

步驟三,對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理，確定區(qū)域中的時(shí)間節(jié)點(diǎn)；對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，確定區(qū)域的空間節(jié)點(diǎn)；

步驟四，確定t＝0時(shí)刻每個(gè)空間節(jié)點(diǎn)的溫度和電場(chǎng)強(qiáng)度的迭代初始值；設(shè)定初始時(shí)刻熱路中的熱端溫度為th(0)，冷端溫度為tl(0)，其中th(0)，tl(0)分別為th(t)和tl(t)在t＝0時(shí)刻的值；

步驟五，對(duì)數(shù)學(xué)模型進(jìn)行離散化處理，建立區(qū)域內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度物理量的代數(shù)方程；步驟六，從溫度和電場(chǎng)強(qiáng)度各空間節(jié)點(diǎn)的迭代初始值出發(fā)，根據(jù)步驟五所建立的節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程，進(jìn)行迭代求解可得出動(dòng)態(tài)溫差發(fā)電系統(tǒng)各空間節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的數(shù)值解，其中，求下一時(shí)間節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的值時(shí)熱端溫度和冷端溫度取上一時(shí)間節(jié)點(diǎn)所在時(shí)刻的值。

步驟三中所述的對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理，確定區(qū)域中的時(shí)間節(jié)點(diǎn)；對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，確定區(qū)域的空間節(jié)點(diǎn)；其特征在于：

(1)對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理，確定區(qū)域中的時(shí)間節(jié)點(diǎn)：

將時(shí)間t均分為(n-1)個(gè)單元，得到n個(gè)時(shí)間節(jié)點(diǎn)，從一個(gè)時(shí)間節(jié)點(diǎn)到下一個(gè)時(shí)間節(jié)點(diǎn)的時(shí)間間隔為δτ稱為時(shí)間步長(zhǎng)，其中第n時(shí)間節(jié)點(diǎn)表示n*δτ時(shí)刻，取n＝0，1，……，n；

(2)對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，確定區(qū)域的空間節(jié)點(diǎn)：

將區(qū)域i均分為(nh-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxi＝li/(nh-1)，其中，li為區(qū)域i的長(zhǎng)度，第i單元格的中心處為空間節(jié)點(diǎn)xi，左邊界和右邊界分別為xi-1/2和xi+1/2，取i＝1，2，……，nh；

將區(qū)域ii均分為(nw-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxii＝lii/(nw-1)，其中，lii為區(qū)域ii的長(zhǎng)度，第j單元的中心處為空間節(jié)點(diǎn)xj，左邊界和右邊界分別為xj-1/2和xj+1/2，取j＝nh+1，nh+2，……，nh+nw；

將區(qū)域iii均分為(nc-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxiii＝liii/(nc-1)，其中l(wèi)iii為區(qū)域iii的長(zhǎng)度，第k單元的中心處為空間節(jié)點(diǎn)xk，左邊界和右邊界分別為xk-1/2和xk+1/2，取k＝nh+nw+1，nh+nw+2，……，nh+nw+nc；

步驟五中所述的對(duì)數(shù)學(xué)模型進(jìn)行離散化處理，建立區(qū)域內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)的溫度和電場(chǎng)強(qiáng)度的代數(shù)方程時(shí)，用空間單元溫度和電場(chǎng)強(qiáng)度的平均值來替代空間節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的值。

本發(fā)明給出一種適用于熱端溫度時(shí)變動(dòng)態(tài)溫差發(fā)電系統(tǒng)的時(shí)域分析方法，該方法給出了一種更普遍適用數(shù)學(xué)計(jì)算模型，首先建立基于熱端溫度時(shí)變動(dòng)態(tài)溫差發(fā)電系統(tǒng)溫差發(fā)電系統(tǒng)的基本物理特性的熱電偏微分方程模型和邊界條件，該模型中所有函數(shù)的假設(shè)與實(shí)際的溫差發(fā)電系統(tǒng)原始的物理特性相聯(lián)系，其次對(duì)系統(tǒng)熱路和電路進(jìn)行了更具體的分析，并進(jìn)行時(shí)間、空間離散化處理，從能量守恒的角度建立區(qū)域內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)物理量的代數(shù)方程，然后從時(shí)間初始值出發(fā)，進(jìn)行迭代求解，最終可得到熱端溫度時(shí)變動(dòng)態(tài)溫差發(fā)電系統(tǒng)內(nèi)部各處溫度和電場(chǎng)強(qiáng)度的數(shù)值解。通過動(dòng)態(tài)溫差發(fā)電內(nèi)部溫度和電場(chǎng)強(qiáng)度的數(shù)值解的分析，可進(jìn)一步對(duì)溫差發(fā)電系統(tǒng)的進(jìn)行更具體的分析和設(shè)計(jì)，故該模型適用于熱端溫度時(shí)變的動(dòng)態(tài)溫差電系統(tǒng)時(shí)域分析，并且結(jié)果更加精確。

附圖說明

圖1為本發(fā)明所述溫差電發(fā)電系統(tǒng)基本結(jié)構(gòu)示意圖；

圖2為本發(fā)明所述溫差發(fā)電系統(tǒng)工作的熱路和電路簡(jiǎn)化示意圖；

圖3為本發(fā)明所述各個(gè)區(qū)域的離散化處理示意圖；

圖4為本發(fā)明的流程圖。

具體實(shí)施方式

為了進(jìn)一步說明本發(fā)明，請(qǐng)參考附圖

如圖1所示，溫差發(fā)電系統(tǒng)由與熱源1接觸的熱端散熱器2，陶瓷片3，導(dǎo)流片4和d對(duì)半導(dǎo)體溫差熱電偶5，與冷源7接觸的冷端散熱器6，陶瓷片3和導(dǎo)流片4組成。為得到較高的電動(dòng)勢(shì)，溫差發(fā)電系統(tǒng)由d對(duì)半導(dǎo)體溫差電偶5經(jīng)電導(dǎo)率較高的導(dǎo)流片4和金屬導(dǎo)線9與負(fù)載電阻8串聯(lián)而成，并夾于導(dǎo)熱性能較好且電絕緣的兩片相互平行的陶瓷片3之間。熱流q流經(jīng)熱端的散熱器2，陶瓷片3和導(dǎo)流片4，d對(duì)半導(dǎo)體溫差熱電偶5，冷端的陶瓷片3，導(dǎo)流片4和冷端散熱器6，從熱源1傳入冷源端7，從熱流通路上來看，所有的溫差電偶為并聯(lián)。而電流j經(jīng)導(dǎo)流片4和金屬導(dǎo)線9流入(出)d對(duì)溫差熱電偶5，與負(fù)載8相連，在電流通路上，所有的熱電偶均為串聯(lián)；熱電耦合只發(fā)出在溫差熱電偶中。

如圖2所示，由于導(dǎo)流片和陶瓷片的導(dǎo)熱性能很好，故忽略其對(duì)熱路的影響，溫差電系統(tǒng)熱路可分成三個(gè)區(qū)域：區(qū)域i，區(qū)域ii和區(qū)域iii，其中區(qū)域i代表熱端散熱器，區(qū)域ii代表d對(duì)半導(dǎo)體溫差熱電偶，區(qū)域iii代表冷端散熱器；其中區(qū)域i和熱源的交界為邊界a，區(qū)域i和區(qū)域ii的交界為邊界b，區(qū)域ii和區(qū)域iii的交界為邊界c，區(qū)域iii和冷源的交界為邊界d。由于導(dǎo)流片和金屬導(dǎo)線的導(dǎo)電性能很好，故忽略其對(duì)電路的影響，系統(tǒng)中電路可分成兩個(gè)區(qū)域：區(qū)域ii和區(qū)域iv；區(qū)域ii代表d對(duì)半導(dǎo)體溫差熱電偶，區(qū)域iv代表負(fù)載電阻；其中區(qū)域ii在熱源端和區(qū)域iv的交界為邊界m，其中區(qū)域ii在冷源端和區(qū)域iv的交界為邊界n。

如圖3所示

(1)時(shí)間的離散化處理：

τ為時(shí)間坐標(biāo)，對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理；

將時(shí)間t均分為(n-1)個(gè)單元，得到n個(gè)時(shí)間節(jié)點(diǎn)，從一個(gè)時(shí)間節(jié)點(diǎn)到下一個(gè)時(shí)間節(jié)點(diǎn)的時(shí)間間隔為δτ稱為時(shí)間步長(zhǎng)，其中第n時(shí)間節(jié)點(diǎn)表示n*δτ時(shí)刻，取n＝0，1，……，n；

(2)一維空間區(qū)域的離散化處理

x為空間坐標(biāo)，對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，由于區(qū)域iv中，物理量的值與空間尺寸無關(guān)，故不對(duì)其進(jìn)行空間離散化處理

將區(qū)域i均分為(nh-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxi＝li/(nh-1)，其中l(wèi)i為區(qū)域i的長(zhǎng)度，第i單元的中心處為空間節(jié)點(diǎn)xi，左邊界和右邊界分別為xi-1/2和xi+1/2，取i＝1，2，……，nh；

將區(qū)域ii，均分為(nw-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxii＝lii/(nw-1)，其中l(wèi)ii為區(qū)域ii的長(zhǎng)度，第j單元的中心處為空間節(jié)點(diǎn)xj，左邊界和右邊界分別為xj-1/2和xj+1/2，取j＝nh+1，nh+2，……，nh+nw；

將區(qū)域iii均分為(nc-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxiii＝liii/(nc-1)，其中l(wèi)iii為區(qū)域iii的長(zhǎng)度，第k單元的中心處為空間節(jié)點(diǎn)xk，左邊界和右邊界分別為xk-1/2和xk+1/2，取k＝nh+nw+1，nh+nw+2，……，nh+nw+nc；

(3)時(shí)間、一維空間離散化后節(jié)點(diǎn)的標(biāo)示

如(n,m)代表了時(shí)間-空間區(qū)域中的位置，n表示時(shí)間節(jié)點(diǎn)位置，m表示空間節(jié)點(diǎn)的位置。m＝i,j,k。

如圖4所示，本發(fā)明的溫差發(fā)電系統(tǒng)動(dòng)態(tài)分析方法包括以下步驟：

步驟01：根據(jù)溫差發(fā)電系統(tǒng)材料的基本物理特性，熱力學(xué)和電學(xué)基本定律，和對(duì)溫差發(fā)電系統(tǒng)電路和熱路的具體分析，導(dǎo)出由d對(duì)半導(dǎo)體溫差熱電偶構(gòu)成的熱端溫度時(shí)變的動(dòng)態(tài)溫差發(fā)電系統(tǒng)熱路和電路的數(shù)學(xué)模型：

(1)由于區(qū)域i和iii只有熱流流過，故根據(jù)熱力學(xué)定律，區(qū)域i和iii的溫度變化偏微分方程的數(shù)學(xué)模型為：

和

其中，ρi和ρiii分別為區(qū)域i和iii物質(zhì)的密度，cvi和cviii分別為區(qū)域i和iii物質(zhì)的定容比熱，ki和kiii分別為區(qū)域i和iii物質(zhì)的熱導(dǎo)率。

(2)熱電耦合發(fā)生在溫差熱電偶中，由于半導(dǎo)體材料中微粒子的作用，熱流密度和電流密度可作為電場(chǎng)強(qiáng)度和溫度梯度的函數(shù)

將其代入電荷存儲(chǔ)方程，高斯定律和能量守恒定律

其中，π＝αt，

其中t為溫度，為溫度梯度；e為電場(chǎng)強(qiáng)度，ρ為物質(zhì)的密度，cv為物質(zhì)的定容比熱，ε為物質(zhì)的介電常數(shù)，j為電流密度，j0為t＝0時(shí)刻的電流密度，α，σ，k，π分別為物質(zhì)的塞貝克系數(shù)，電導(dǎo)率，熱導(dǎo)率和珀?duì)柼禂?shù)。

整理可得，區(qū)域ii中有：

其中t為溫度，e為電場(chǎng)強(qiáng)度，ρii為區(qū)域ii物質(zhì)的密度，cvii為區(qū)域ii物質(zhì)的定容比熱，εii為區(qū)域ii物質(zhì)的介電常數(shù)，j為電流密度，j0為t＝0時(shí)刻的電流密度，αii，σii，kii，為區(qū)域ii物質(zhì)的總等效塞貝克系數(shù)，總等效電導(dǎo)率，總等效熱導(dǎo)率。由于區(qū)域ii由d對(duì)p型和n型半導(dǎo)體溫差熱電偶在電學(xué)上串聯(lián)，在熱學(xué)上并聯(lián)而成，故αii＝d*(αp-αn)，kii＝m*(kp+kn)；αp，σp，kp，分別為溫差熱電偶中p型半導(dǎo)體材料的塞貝克系數(shù)，電導(dǎo)率和熱導(dǎo)率，αn，σn，kn，分別為溫差熱電偶中n型半導(dǎo)體材料的塞貝克系數(shù)，電導(dǎo)率和熱導(dǎo)率。

(3)由于區(qū)域iv中負(fù)載電阻rl只有電流流過，根據(jù)電路基本定律有：

區(qū)域iv中，

σiv為負(fù)載電阻的電導(dǎo)率。

(4)整理可得，熱路各個(gè)區(qū)域控制方程為：

區(qū)域i中，

區(qū)域ii中，

區(qū)域iii中，

電路各個(gè)區(qū)域的控制方程為：

區(qū)域ii中，

區(qū)域iv中，

步驟02：根據(jù)實(shí)際工況條件確定熱端溫度時(shí)變動(dòng)態(tài)溫差發(fā)電系統(tǒng)的熱路和電路的邊界條件：

(1)熱路中的邊界條件為：

邊界a:根據(jù)熱源提供能量的形式，熱端等效為溫度為th(t)；

邊界b:通過邊界處的熱流密度為連續(xù)的，溫度為連續(xù)的；

邊界c:通過邊界處的熱流密度為連續(xù)的，溫度為連續(xù)的；

邊界d：根據(jù)冷源端散熱的形式，冷端等效為溫度為tl(t)；

(2)電路中的邊界條件為：

邊界m:通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的；

邊界n:通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的；

步驟03：對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理，確定區(qū)域中的時(shí)間節(jié)點(diǎn)，對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，確定區(qū)域中的空間節(jié)點(diǎn)；

(1)時(shí)間離散化處理

τ為時(shí)間坐標(biāo)，對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii，區(qū)域iv進(jìn)行時(shí)間離散化處理

將時(shí)間t均分為(n-1)個(gè)單元，得到n個(gè)時(shí)間層，從一個(gè)時(shí)間層到下一個(gè)時(shí)間層的時(shí)間間隔為δτ稱為時(shí)間步長(zhǎng)，其中第n時(shí)間層表示n*δτ時(shí)刻，取n＝0，1，……，n；

(2)一維空間區(qū)域的離散化處理

x為空間坐標(biāo)，對(duì)區(qū)域i，區(qū)域ii，區(qū)域iii進(jìn)行一維空間離散化處理，由于區(qū)域iv中，物理量的值與空間尺寸無關(guān)，故無需對(duì)其進(jìn)行空間離散化處理

將區(qū)域i均分為(nh-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxi＝li/(nh-1)，其中l(wèi)i為區(qū)域i的長(zhǎng)度，第i單元格的中心處為空間節(jié)點(diǎn)xi，左邊界和右邊界分別為xi-1/2和xi+1/2，取i＝1，2，……，nh；

將區(qū)域ii均分為(nw-1)個(gè)單元,則每個(gè)單元的長(zhǎng)度為δxii＝lii/(nw-1)，其中l(wèi)ii為區(qū)域ii的長(zhǎng)度，第j單元格的中心處為空間節(jié)點(diǎn)xj，左邊界和右邊界分別為xj-1/2和xj+1/2，取j＝nh+1，nh+2，……，nh+nw；

將區(qū)域iii均分為(nc-1)個(gè)單元,則每個(gè)單元格的長(zhǎng)度為δxiii＝liii/(nc-1)，其中l(wèi)iii為區(qū)域iii的長(zhǎng)度，第k單元格的中心處為空間節(jié)點(diǎn)xk，左邊界和右邊界分別為xk-1/2和xk+1/2，取k＝nh+nw+1，nh+nw+2，……，nh+nw+nc；

(3)時(shí)間、一維空間離散后物理量的標(biāo)示

eⁿ表示物理量e在第n時(shí)間層的值；tm表示物理量t在m空間節(jié)點(diǎn)的值；m＝i,j,k。如(n,m)代表了時(shí)間-空間區(qū)域中的位置，n表示時(shí)間節(jié)點(diǎn)位置，m表示空間節(jié)點(diǎn)的位置。表示物理量t第n時(shí)間m空間節(jié)點(diǎn)的值，表示在n*δτ時(shí)刻第m單元溫度平均值，m＝i,j,k；同樣地，該節(jié)點(diǎn)的電場(chǎng)強(qiáng)度記為表示在n*δτ時(shí)刻第m單元電場(chǎng)強(qiáng)度的平均值，m＝j(luò)。步驟04：確定t＝0時(shí)刻每個(gè)節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的迭代初始值；

首先確定溫度t和電場(chǎng)強(qiáng)度e在t＝0時(shí)刻每個(gè)空間節(jié)點(diǎn)的初始值和

設(shè)定初始時(shí)刻熱路中的熱端溫度為th0＝th(0)，冷端溫度為tl0＝tl(0)，

則初始時(shí)刻熱流密度相同設(shè)定為q0＝(th(0)-tl(0))/(1/ki+1/kii+1/kiii)，根據(jù)熱路方程可得初始時(shí)刻各空間節(jié)點(diǎn)的溫度初始值為：

邊界處溫度初始值為：

根據(jù)電路方程可得到電流密度的初始值為

進(jìn)而可得到初始時(shí)刻空間內(nèi)部節(jié)點(diǎn)和邊界處電場(chǎng)強(qiáng)度的初始值

其中，由于區(qū)域ii在熱路的邊界c和d與區(qū)域ii在電路的邊界m，n重合。故有

步驟05：對(duì)數(shù)學(xué)模型進(jìn)行離散化處理，建立區(qū)域內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程；(1)建立區(qū)域內(nèi)部節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程

用空間單元物理量的平均值來代替節(jié)點(diǎn)物理量的值，根據(jù)泰勒級(jí)數(shù)展開式，有：

對(duì)區(qū)域i，對(duì)偏微分方程進(jìn)行一維離散化處理可得

則：

同樣地，

對(duì)區(qū)域iii，對(duì)偏微分方程進(jìn)行一維離散化處理可得

對(duì)區(qū)域ii，對(duì)偏微分方程一維形式為

離散化處理可得

則

對(duì)時(shí)間節(jié)點(diǎn)有

最終，可整理出區(qū)域內(nèi)部節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程為：

區(qū)域i中

區(qū)域ii

區(qū)域iii

區(qū)域iv中，

(2)建立區(qū)域邊界節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程：

在邊界a處,根據(jù)熱源提供能量的形式，熱端等效為溫度為th(t)，則在n*δτ時(shí)刻的熱端溫度設(shè)為th(n*δτ)，這里，取n＝1，……，n；

在邊界d處，根據(jù)冷端散熱方式的形式，冷端等效為溫度為tl(t)，則在n*δτ時(shí)刻的冷端溫度為tl(n*δτ)，這里，取n＝1，……，n；

對(duì)于邊界b，根據(jù)熱路的邊界條件任一時(shí)刻熱流密度是連續(xù)的，任一時(shí)刻溫度是連續(xù)的

有

故有

離散化處理可得，這里i＝nh，j＝nh+1

對(duì)于邊界m，根據(jù)電路方程和電路中邊界條件任一時(shí)刻通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的，有

故有

由于在區(qū)域ii中電路的邊界m和熱路的邊界b為同一位置，

故其中

聯(lián)立方程(1-1)和(1-2)可得邊界b處

同樣地，對(duì)于邊界c，根據(jù)熱路的邊界條件任一時(shí)刻熱流密度是連續(xù)的，任一時(shí)刻溫度是連續(xù)的

有

故有

離散化處理可得，這里j＝nh+nw，k＝nh+nw+1

對(duì)于邊界n，根據(jù)電路方程和電路中邊界條件任一時(shí)刻通過邊界處的電流密度為連續(xù)的，電場(chǎng)強(qiáng)度是連續(xù)的。

有

故有

由于在區(qū)域ii中電路的邊界n和熱路的邊界c為同一位置，

故其中

聯(lián)立方程(2-1)和(2-2)可得邊界c處

最終，可整理出區(qū)域各節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的離散化代數(shù)方程為：

區(qū)域i中

區(qū)域ii中，

其中，

(i＝nh，j＝nh+1)

(j＝nh+nw，k＝nh+nw+1)

步驟06，從溫度和電場(chǎng)強(qiáng)度各空間節(jié)點(diǎn)的迭代初始值出發(fā)，根據(jù)步驟05所建立的節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的代數(shù)方程，進(jìn)行迭代求解可得出動(dòng)態(tài)溫差發(fā)電系統(tǒng)各節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的數(shù)值解，其中，求(n+1)時(shí)間節(jié)點(diǎn)溫度和電場(chǎng)強(qiáng)度的數(shù)值解時(shí)熱冷端的邊界條件取n*δτ時(shí)刻的值。