技术领域
[0001] 本
发明涉及
半导体温差发电技术领域，尤其涉及一种温差热电偶材料分段动态温差发电系统时域分析方法。
背景技术
[0002] 1821年德国科学家Seebeck发现，在两种不同的金属构成的闭合回路中，当两个接头存在温差时，回路将产生
电流，这个现象被称为
塞贝克效应。温差发电就是利用塞贝克效应，通过在热电材料两端维持一定的温差，从而产生一定的
电压和电功率输出。通过研究发现，半导体材料的热电优值较大，目前人们在温差系统中应用的热电材料都是半导体材料，所以又称为半导体温差发电。
[0003] 随着人们对
能源危机的认识，温差发电技术可以利用自然界中存在的大量的温差以及工业余热，具有良好的综合社会经济效益。同时，随着人们对空间探索兴趣，医用物理学的进展，和大规模无线
传感器的应用，需要开发一类能够自身供给
能量并无需照看的电源系统，显然，温差发电对这些应用极为合适。温差发电作为一种全固态能量转换方式，具有无介质泄露，无噪声，性能可靠，维护少等优点的绿色环保能源，在微型能源、低品位能源、废能源利用方面的应用价值越来越明显。因而，尽快实现温差发电技术及其应用的产业化具有重要的现实意义。而温差发电系统的分析和设计依赖于实际的工作状况，故针对不同的实际情况进行合理的温差发电系统模型是对温差电系统的各种参数进行更具体的分析和设计的
基础。
[0004] 中国
专利文献CN201410189306.6公开了一种温差发电系统电动势计算方法，在该专利中采用的对温差电系统进行分析计算的数学模型，对实际温差发电系统进行了如下的近似：
[0005] 认为温差电系统中的物理特性参数是恒定值，即为常数；
[0006] 认为边界处的热流或
温度是时不变的，即为常数，
[0007] 故专利文献CN201410189306.6中温差电系统电动势的计算方法，将不同物质交界处的边界条件，简单近似为连接处的
流体平衡，只适用于温差发电系统是处于常物性的，边界处可简单近似为平衡态情况下，稳态线性温差电系统的分析，而对于温差发电系统处于非平衡态时是不适用的，即不能用于分析动态温差发电系统的瞬态过程。
[0008] 而有时候实际的温差电系统中，许多材料的最大优值温度范围通常较窄，很难使材料达到最佳掺杂分布，也难于
覆盖整个应用温区，为了提高整个温差电偶在整个工作温区内的平均优值，采用分段结构以提高器件的发电性能，当温差较小时，其中的各段可采用同一种材料，进行不同的掺杂而获得，但对于温差较大的情形，由于材料最大优值范围的熔点等限制，其中各段需要采用不同的材料构成。因此，实际的温差电偶是由几种材料组成时，在不同的分段区域物质的物理特性参数很显然是不相同的，不再是恒定值。
发明内容
[0009] 本发明的目的是提供一种温差热电偶材料分段动态温差发电系统时域分析方法，通过本发明中的方法得到动态温差发电内部温度和
电场强度的数值解，可进一步对热端热流时变动态温差发电系统的进行更具体的分析和设计，并且结果更加精确。
[0010] 本发明采用的技术方案：
[0011] 一种温差热电偶分段动态温差发电系统时域分析方法，具体包括以下几个步骤：
[0012] 步骤一，根据温差发电系统材料的基本物理特性，热
力学和电学基本定律，导出由d对分b段半导体温差电偶构成的温差热电偶分段的动态温差发电系统热路和
电路的数学模型：
[0013] 热路中各个区域的数学模型为
[0014] 区域I中，
[0015] 区域IIu中，
[0016] 区域III中，
[0017] 电路各个区域的数学模型为：
[0018] 区域IIu中，
[0019]
[0020] 区域IV中，
[0021] 其中，温差发电系统热路分成三个区域：区域I，区域II和区域III，区域I代表热端
散热器，区域II代表d对材料分b段的半导体温差热电偶，区域III代表冷端
散热器；温差发电系统中电路分成两个区域：区域II和区域IV；区域II代表d对材料分b段的半导体温差热电偶，区域IV代表负载
电阻RL；区域II中b段不同材料分别对应区域II1，II2,…IIu，即区域IIu对应区域II中第u段材料，u＝1，2，……，b；
[0022] 其中，ρI和ρIII分别为区域I和III物质的
密度，CvI和CvIII分别为区域I和III物质的定容
比热，kI和kIII分别为区域I和III物质的热导率。其中T为温度， 为温度梯度，E为电场强度；ρIIu分别为区域IIu中物质的密度，CvIIu为区域IIu中物质的定容比热，εIIu为区域IIu中物质的
介电常数，J为电路中的电流密度，电流J0为t＝0时刻通过电路的电流密度，αIIu，σIIu，kIIu，分别为区域IIu中物质的总等效
塞贝克系数，总等效电导率，总等效热导率；其中αIIu＝d*(αPu-αNu)， kII＝d*(kPu+kNu)；αPu，σPu，kPu，分别为区域IIu中P型半导体材料的塞贝克系数，电导率和热导率，αNu，σNu，kNu，分别为区域IIu中N型半导体材料的塞贝克系数，电导率和热导率；σIV为负载电阻的电导率。
[0023] 步骤二，根据实际工况条件确定热路和电路的边界条件：
[0024] 热路中的边界条件为：
[0025] 边界A:根据热源提供能量的形式，热端等效为热流密度为q0；
[0026] 边界B:通过边界处的热流密度为连续的，温度为连续的；
[0027] 边界Pu:通过边界处的热流密度为连续的，温度为连续的；
[0028] 边界C:通过边界处的热流密度为连续的，温度为连续的；
[0029] 边界D：根据冷源端散热的形式，冷端等效为温度为Tl0；
[0030] 电路中的边界条件为：
[0031] 边界M:通过边界处的电流密度为连续的，电场强度是连续的；
[0032] 边界Pu:通过边界处的电流密度为连续的，电场强度为连续的；
[0033] 边界N:通过边界处的电流密度为连续的，电场强度是连续的；
[0034] 其中热路中区域I和热源的交界为边界A，区域I和区域II的交界为边界B，区域IIu和II(u+1)的交界为边界Pu；区域II和区域III的交界为边界C，区域III和冷源的交界为边界D；其中电路中区域II在热源端和区域IV的交界为边界M，区域IIu和II(u+1)的交界为边界Pu；区域II在冷源端和区域IV的交界为边界N；取u＝1，2，……，b-1；
[0035] 步骤三,对区域I，区域II，区域III，区域IV进行时间离散化处理，确定区域中的时间
节点；对区域I，区域II，区域III进行一维空间离散化处理，确定区域的空间节点；
[0036] 步骤四，确定t＝0时刻每个空间节点的温度和电场强度的
迭代初始值；
[0037] 步骤五，对数学模型进行离散化处理，建立区域内部节点和边界节点温度和电场强度的代数方程；
[0038] 步骤六，从温度和电场强度各空间节点的迭代初始值出发，根据步骤五所建立的节点温度和电场强度的代数方程，进行迭代求解可得出温差热电偶材料分段的动态温差发电系统各节点温度和电场强度的数值解；
[0039] 步骤三中所述的对区域I，区域II，区域III，区域IV进行时间离散化处理，确定区域中的时间节点；对区域I，区域II，区域III进行一维空间离散化处理，确定区域的空间节点；其特征在于：
[0040] (1)对区域I，区域II，区域III，区域IV进行时间离散化处理，确定区域中的时间节点：
[0041] 将时间t均分为(N-1)个单元，得到N个时间节点，从一个时间节点到下一个时间节点的时间间隔为Δτ称为时间步长，其中 第n时间节点表示n*Δτ时刻，取n＝0，1，……，N；
[0042] (2)对区域I，区域II，区域III进行一维空间离散化处理，确定区域的空间节点：
[0043] 将区域I均分为(Nh-1)个单元,则每个单元格的长度为ΔxI＝LI/(Nh-1)，其中，LI为区域I的长度，得到Nh个空间节点，第i单元格的中心处为空间节点xi，左边界和右边界分别为xi-1/2和xi+1/2，取i＝1，2，……，Nh；
[0044] 区域II中对应的b段不同材料分别为区域II1，II2,…IIu，即区域IIu对应区域II中第u段材料，各段区域分别进行空间离散化处理，将区域IIu均分为(Nu-1)个单元,则区域IIu中每个单元格的长度为ΔxIIu＝LIIu/(Nu-1)，其中LIIu为区域IIu的长度，第ju单元格的中心处为空间节点xju，左边界和右边界分别为xju-1/2和xju+1/2，取ju＝Nh+N1+N(u-1)+1，……，Nh+N1+N(u-1)+Nu，取u＝1，2，……，b；令Nw＝N1+N2+...+Nb；
[0045] 将区域III均分为(Nc-1)个单元,则每个单元格的长度为ΔxIII＝LIII/(Nc-1)，其中LIII为区域III的长度，第k单元格的中心处为空间节点xk，左边界和右边界分别为xk-1/2和xk+1/2，取k＝Nh+Nw+1，Nh+Nw+2，……，Nh+Nw+Nc；
[0046] 步骤五中所述的对数学模型进行离散化处理，建立区域内部节点和边界节点的温度和电场强度的代数方程时，用空间单元温度和电场强度的平均值来替代空间节点温度和电场强度的值。
[0047] 本发明给出一种适用于温差热电偶材料分段的动态温差发电系统的时域分析方法，该方法给出了一种更普遍适用数学计算模型，首先建立基于温差热电偶材料分段的动态温差发电系统的基本物理特性的热电偏微分方程模型和边界条件，该模型中所有函数的假设与实际的温差发电系统原始的物理特性相联系，其次对系统热路和电路进行了更具体的分析，并进行时间、空间离散化处理，从能量守恒的
角度建立区域内部节点和边界节点物理量的代数方程，然后从时间初始值出发，进行迭代求解，最终可得到温差热电偶材料分段的动态温差发电系统内部各处温度和电场强度的数值解。通过动态温差发电内部温度和电场强度的数值解的分析，可进一步对温差发电系统的进行更具体的分析和设计，故该模型适用于温差热电偶材料分段的动态温差发电系统时域分析，并且结果更加精确。
附图说明
[0048] 图1为本发明所述温差电发电系统基本结构示意图；
[0049] 图2为本发明所述温差发电系统工作的热路和电路简化示意图；
[0050] 图3为本发明所述各个区域的离散化处理示意图；
[0051] 图4为本发明所述区域II中材料分段的温差热电偶离散化处理示意图；
具体实施方式
[0053] 为了进一步说明本发明，请参考附图
[0054] 如图1所示，温差发电系统由与热源1
接触的热端散热器2，陶瓷片3，导流片4和d对半导体温差热电偶5，与冷源7接触的冷端散热器6，陶瓷片3和导流片4组成。为得到较高的电动势，温差发电系统由d对半导体温差电偶5经电导率较高的导流片4和金属
导线9与负载电阻8
串联而成，并夹于导热性能较好且电绝缘的两片相互平行的陶瓷片3之间。热流q流经热端的散热器2，陶瓷片3和导流片4，d对半导体温差热电偶5，冷端的陶瓷片3，导流片4和冷端散热器6，从热源1传入冷源端7，从热流通路上来看，所有的温差电偶为并联。而电流J经导流片4和金属导线9流入(出)d对温差热电偶5，与负载8相连，在电流通路上，所有的热电偶均为串联；故热电耦合只发出在温差热电偶中。
[0055] 图2为温差发电系统工作的热路和电路简化示意图，由于导流片和陶瓷片的导热性能很好，故忽略其对热路的影响，温差电系统热路可分成三个区域：区域I，区域II和区域III，其中区域I代表热端散热器，区域II代表d对半导体温差热电偶，区域III代表冷端散热器；其中区域I和热源的交界为边界A，区域I和区域II的交界为边界B，区域II和区域III的交界为边界C，区域III和冷源的交界为边界D。由于导流片和金属导线的
导电性能很好，故忽略其对电路的影响，系统中电路可分成两个区域：区域II和区域IV；区域II代表d对半导体温差热电偶，区域IV代表负载电阻；其中区域II在热源端和区域IV的交界为边界M，其中区域II在冷源端和区域IV的交界为边界N。
[0056] 如图3和图4所示，各个区域的离散化处理如下：
[0057] (1)时间离散化处理
[0058] τ为时间坐标，对区域I，区域II，区域III，区域IV进行时间离散化处理[0059] 将时间t均分为(N-1)个单元，得到N个时间层，从一个时间层到下一个时间层的时间间隔为Δτ称为时间步长，其中 第n时间层表示n*Δτ时刻，取n＝0，1，……，N；En表示物理量E在第n时间层的值；
[0060] (2)一维空间区域的离散化处理
[0061] x为空间坐标，对区域I，区域II，区域III进行一维空间离散化处理，由于区域IV中，物理量的值与空间尺寸无关，故不对其进行空间离散化处理
[0062] 将区域I均分为(Nh-1)个单元,则每个单元格的长度为ΔxI＝LI/(Nh-1)，其中LI为区域I的长度，得到Nh个空间节点，第i单元格的中心处为空间节点xi，左边界和右边界分别为xi-1/2和xi+1/2，取i＝1，2，……，Nh；
[0063] 区域II中对应的b段不同材料分别为区域II1，II2,…IIu，即区域IIu对应区域II中第u段材料，各段区域分别进行空间离散化处理，将区域IIu均分为(Nu-1)个单元,则区域IIu中每个单元格的长度为ΔxIIu＝LIIu/(Nu-1)，其中LIIu为区域IIu的长度，第ju单元格的中心处为空间节点xju，左边界和右边界分别为xju-1/2和xju+1/2，取ju＝Nh+N1+N(u-1)+1，……，Nh+N1+N(u-1)+Nu，取u＝1，2，……，b；令Nw＝N1+N2+...+Nb；
[0064] 将区域III均分为(Nc-1)个单元,则每个单元格的长度为ΔxIII＝LIII/(Nc-1)，其中LIII为区域III的长度，第k单元格的中心处为空间节点xk，左边界和右边界分别为xk-1/2和xk+1/2，取k＝Nh+Nw+1，Nh+Nw+2，……，Nh+Nw+Nc；Tm表示物理量T在m空间节点的值；m＝i,ju,k取u＝1，2，……，b。
[0065] (3)时间、一维空间的离散化后节点物理量的标示
[0066] 如(n,m)代表了时间-空间区域中的
位置，n表示时间节点位置，m表示空间节点的位置。m＝i,ju,k；取u＝1，2，……，b。
[0067] 如图5所示，本发明的温差热电偶材料分段的动态温差发电系统包括以下步骤：
[0068] 步骤01：根据温差热电偶材料分段的动态温差发电系统的基本物理特性，
热力学和电学基本定律，和对温差发电系统电路和热路的具体分析，导出由d对半导体温差热电偶构成的温差发电系统热路和电路的数学模型；
[0069] (1)由于区域I和III只有热流流过，故根据热力学定律，区域I和III的温度变化偏微分方程的数学模型为：
[0070]
[0071] 其中，ρI和ρIII分别为区域I和III物质的密度，CvI和CvIII分别为区域I和III物质的定容比热，kI和kIII分别为区域I和III物质的热导率。
[0072] (2)热电耦合发生在温差热电偶中，由于半导体材料中微粒子的作用，热流密度和电流密度可作为电场强度和温度梯度的函数
[0073]
[0074]
[0075] 将其代入电荷存储方程，高斯定律和能量守恒定律
[0076]
[0077]
[0078] 其中，π＝αT，
[0079] 其中T为温度， 为温度梯度，E为电场强度，ρ分别为物质密度，Cv为物质的定容比热，ε为物质的介电常数，J为电流密度，α，σ，k，π为物质的塞贝克系数，电导率，热导率和珀尔帖系数。
[0080] 区域II代表d对材料分b段的半导体温差热电偶，区域II中b段不同材料分别对应区域II1，II2,…IIu，即区域IIu对应区域II中第u段材料，u＝1，2，……，b；故
整理可得[0081] 区域II中有：
[0082]
[0083]
[0084]
[0085] ρIIu分别为区域IIu中物质的密度，CvIIu为区域IIu中物质的定容比热，εIIu为区域IIu中物质的介电常数，J为电路中的电流密度，电流J0为t＝0时刻通过电路的电流密度，αIIu，σIIu，kIIu，为区域IIu中物质的总等效塞贝克系数，总等效电导率，总等效热导率；其中αIIu＝d*(αPu-αNu)， kII＝d*(kPu+kNu)；αPu，σPu，kPu，分别为区域IIu中P型半导体材料的塞贝克系数，电导率和热导率，αNu，σNu，kNu，分别为区域IIu中N型半导体材料的塞贝克系数，电导率和热导率；σIV为负载电阻的电导率。
[0086] (3)由于区域IV中负载电阻RL只有电流流过，根据电路基本定律有：
[0087] 区域IV中，
[0088] (4)故整理可得，热路中各个区域的数学模型为
[0089] 区域I中，
[0090] 区域IIu中，
[0091] 区域III中，
[0092] 电路各个区域的数学模型为：
[0093] 区域IIu中，
[0094]
[0095] 区域IV中，
[0096] u＝1，2，……，b；
[0097] 步骤02：根据实际工况条件确定温差热电偶材料分段的动态温差发电系统热路和电路的边界条件；
[0098] 热路中的边界条件为：
[0099] 边界A:根据热源提供能量的形式，热端等效为热流密度为q0；
[0100] 边界B:通过边界处的热流密度为连续的，温度为连续的；
[0101] 边界Pu:通过边界处的热流密度为连续的，温度为连续的；
[0102] 边界C:通过边界处的热流密度为连续的，温度为连续的；
[0103] 边界D：根据冷源端散热的形式，冷端等效为温度为Tl0；
[0104] 电路中的边界条件为：
[0105] 边界M:通过边界处的电流密度为连续的，电场强度是连续的；
[0106] 边界Pu:通过边界处的电流密度为连续的，电场强度为连续的；
[0107] 边界N:通过边界处的电流密度为连续的，电场强度是连续的；
[0108] 其中热路中区域I和热源的交界为边界A，区域I和区域II的交界为边界B，区域IIu和II(u+1)的交界为边界Pu；区域II和区域III的交界为边界C，区域III和冷源的交界为边界D；其中电路中区域II在热源端和区域IV的交界为边界M，区域IIu和II(u+1)的交界为边界Pu；区域II在冷源端和区域IV的交界为边界N；取u＝1，2，……，b-1；
[0109] 步骤03：对区域I，区域II，区域III，区域IV进行时间离散化处理，确定区域中的时间节点，对区域I，区域II，区域III进行一维空间离散化处理，确定区域中的空间节点；
[0110] (1)时间离散化处理
[0111] τ为时间坐标，对区域I，区域II，区域III，区域IV进行时间离散化处理[0112] 将时间t分割为(N-1)个单元，得到N个时间节点，从一个时间节点到下一个时间节点的时间间隔为Δτ称为时间步长，其中 第n时间节点表示n*Δτ时刻，取n＝0，n1，……，N；E表示物理量E在第n时间层的值；
[0113] (2)一维空间区域的离散化处理
[0114] x为空间坐标，对区域I，区域II，区域III进行一维空间离散化处理，由于区域IV中，物理量的值与空间尺寸无关，故不对其进行空间离散化处理
[0115] 将区域I分割为(Nh-1)个单元,则每个单元格的长度为ΔxI＝LI/(Nh-1)，其中LI为区域I的长度，得到Nh个空间节点，第i单元格的中心处为空间节点xi，左边界和右边界分别为xi-1/2和xi+1/2，取i＝1，2，……，Nh；
[0116] 将区域II对应b段不同材料分为b部分，即区域IIu对应区域II中第u段材料，各段区域分别进行空间离散化处理，将区域IIu均分为(Nu-1)个单元,则每个单元格的长度为ΔxIIu＝LIIu/(Nu-1)，其中LIIu为区域IIu的长度，第ju单元格的中心处为空间节点xju，左边界和右边界分别为xju-1/2和xju+1/2，取ju＝Nh+N1+N(u-1)+1，……，Nh+N1+N(u-1)+Nu，取u＝1，2，……，b；令Nw＝N1+N2+...+Nb；
[0117] 将区域III分割为(Nc-1)个单元,则每个单元格的长度为ΔxIII＝LIII/(Nc-1)，其中LIII为区域III的长度，第k单元格的中心处为空间节点xk，左边界和右边界分别为xk-1/2和xk+1/2，取k＝Nh+Nw+1，Nh+Nw+2，……，Nh+Nw+Nc；Tm表示物理量T在m空间节点的值；m＝i,ju,k,取u＝1，2，……，b。
[0118] (3)时间、一维空间离散化后物理量的标示
[0119] 如(n,m)代表了时间-空间区域中的位置，n表示时间节点位置，m表示空间节点的位置。Tmn表示第n时间层m空间节点的温度，， 表示在n*Δτ时刻第m单元温度平均值，m＝i,ju,k；同样地，该节点的电场强度记为 表示在n*Δτ
时刻第m单元电场强度的平均值，m＝ju,取u＝1，2，……，b。
[0120] 步骤04：确定t＝0时刻每个节点温度和电场强度的迭代初始值；
[0121] 首先确定温度T和电场强度E在t＝0时刻每个单元的初始值
[0122] 根据热路方程可得初始时刻各单元的温度初始值为：
[0123]
[0124]
[0125]
[0126]
[0127] 边界处初始值为
[0128]
[0129]
[0130]
[0131]
[0132]
[0133] 根据电路方程可得到电流密度的初始值为
[0134]
[0135] 进而可得到初始时刻单元的电场强度的初始值
[0136]
[0137]
[0138]
[0139] 其中，由于区域II在热路的边界C和D与区域II在电路的边界M，N重合。故有[0140] 步骤05：对数学模型进行离散化处理，建立区域内部节点和边界节点温度和电场强度的代数方程；
[0141] (1)建立区域内部节点温度和电场强度的代数方程
[0142] 用空间单元物理量的平均值来代替节点物理量的值，根据泰勒级数展开式，有：
[0143]
[0144]
[0145]
[0146]
[0147]
[0148] 对区域I，对偏微分方程进行一维离散化处理可得
[0149]
[0150]
[0151]
[0152] 则：
[0153]
[0154] 同样地，
[0155] 对区域III，对偏微分方程进行一维离散化处理可得
[0156]
[0157] 对区域IIu，对偏微分方程一维形式为
[0158]
[0159]
[0160] 离散化处理可得
[0161]
[0162]
[0163]
[0164]
[0165] 则
[0166]
[0167]
[0168]
[0169] 对时间节点有
[0170]
[0171]
[0172] 最终，可整理出区域内部节点温度和电场强度的代数方程为：
[0173] 区域I中
[0174]
[0175] 区域IIu中，
[0176]
[0177]
[0178]
[0179] 区域III
[0180]
[0181] 区域IV中，
[0182] (2)建立区域边界节点物理量的代数方程：
[0183] 在边界A处,根据热源提供能量的形式，热端等效为热流密度为q0，
[0184]
[0185] 在边界D处，根据冷端散热方式的形式，冷端等效为温度为Tl0，
[0186]
[0187] 边界B，区域I与区域II的交界，即为区域I与区域II1的交界，根据热路的边界条件任一时刻热流密度是连续的，任一时刻温度是连续的
[0188] 有
[0189] 故有
[0190]
[0191] 离散化处理可得，这里i＝Nh，j1＝Nh+1
[0192]
[0193]
[0194] 对于边界M，即为区域IV与区域II1的交界，根据电路方程和电路中边界条件任一时刻通过边界处的电流密度为连续的，电场强度是连续的。
[0195] 有
[0196] 故有
[0197] 由于在区域II中电路的边界M和热路的边界B为同一位置，
[0198] 故其中
[0199]
[0200] 联立方程(1-1)和(1-2)可得边界B处
[0201]
[0202]
[0203] 同样地，对于边界C，区域III与区域II的交界，即为区域III与区域IIb的交界根据热路的边界条件任一时刻热流密度是连续的，任一时刻温度是连续的
[0204] 有
[0205] 故有
[0206]
[0207] 离散化处理可得，这里jb＝Nh+Nw，k＝Nh+Nw+1
[0208]
[0209]
[0210] 对于边界N，即为区域IV与区域IIb的交界根据电路方程和电路中边界条件任一时刻通过边界处的电流密度为连续的，电场强度是连续的。
[0211] 有
[0212] 故有
[0213] 由于在区域II中电路的边界N和热路的边界C为同一位置，
[0214] 故其中
[0215]
[0216] 联立方程(2-1)和(2-2)可得边界C处
[0217]
[0218]
[0219] 同样地，对于边界Pu，有
[0220] 根据热路的控制方程和边界条件任一时刻热流密度是连续的，任一时刻温度是连续的，
[0221]
[0222] 根据电路的控制方程和边界条件任一时刻通过边界处的电流密度为连续的，电场强度是连续的。
[0223]
[0224] 故有
[0225]
[0226]
[0227] 离散化处理可得，这里ju＝Nh+N1+...+Nu，j(u+1)＝Nh+N1+...+Nu+1，取u＝1，2，……，b；
[0228]
[0229]
[0230]
[0231]
[0232] 联立方程(3-1)和(3-2)可得边界Pu处
[0233]
[0234]
[0235] 最终，可整理出区域各节点温度和电场强度的离散化代数方程为：
[0236] 区域I中
[0237]
[0238]
[0239]
[0240] 区域II中，
[0241]
[0242]
[0243]
[0244]
[0245]
[0246]
[0247]
[0248]
[0249] 其中，
[0250]
[0251]
[0252]
[0253]
[0254]
[0255]
[0256] 步骤06，从温度和电场强度各空间节点的迭代初始值出发，根据步骤05所建立的节点温度和电场强度的代数方程，进行迭代求解可得出动态温差发电系统各空间节点温度和电场强度的数值解。
