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