基于TPS法测定乙醇热导率：Python绘图计算问题求助

我正在用TPS法（热盘传感器）测定乙醇热导率，需要用Python绘制无量纲时间τ的函数∆D(τ)，以及∆T-D(τ)曲线。根据参考文献《Rapid thermal conductivity measurement with a hot disk sensor. Part 1. Theoretical considerations》，当τ取值正确时，∆T-D(τ)应为直线，以此计算热导率。但编写代码后，无论热扩散率κ取何值，该曲线始终为直线，怀疑修正贝塞尔函数的使用存在问题。

实验参数 #

• 时间范围：0.015s至15s，共1000个点
• κ取值范围：0.00000001至0.3，共1000个点

问题代码 #

from numpy import *
import scipy.special
from scipy.integrate import quad
from matplotlib.pyplot import *
def integer1(sigma):
    return 1/(sigma**2)
 
tini = 0.015
tfin = 15
time = linspace(tini,tfin,num=1000)
n=7 # number of concentric circles of sensor 
L = 1
L0 = np.i0(l*k)/(2*thau**2*n**2)
P0 = 0.1 #power
k = 1 #thermal diffusivity
a = 0.000958 # radius of biggest ring 
λ = 0.169 #thermal conductivity of ethanol (im not sure if this is ok)
x=linspace(0.00000001,0.3,1000)
for K in range (0,len(x)):
    # print (x[K])
    
    theta = a**2/x[K]
    Tlist = []
    Dlist = []
    
    for t in time:
    
        thau = sqrt(t/theta)
        
        som = 0
        for l in range(L,n):
            for k in range(1,n):
                
                som += l*k*exp((-l**2+-k**2)/(4*thau**2*n**2))*L0
        
        I = quad(integer1, 0, thau)
        
        D = ((n*(n+1))**-2)*I[0]*som
        
        T = (P0/(pi**(3/2)*a*λ))*D
        
        Tlist.append(T)
        Dlist.append(D)
        
    figure(1)
    plot(Dlist,Tlist)
show()

核心问题排查 #

1. 强制线性导致曲线恒直
   代码中T = (P0/(pi**(3/2)*a*λ))*D直接让∆T与D(τ)成固定线性比例，这是人为制造的线性关系，和κ取值无关——无论D怎么变，T都是D的固定倍数，必然导致曲线始终为直线。这是最核心的逻辑错误，并非贝塞尔函数的问题。

2. 修正贝塞尔函数调用错误

   • 全局定义的L0 = np.i0(l*k)/(2*thau**2*n**2)存在未定义变量问题：l和thau此时还未进入循环，属于提前引用；
   • 全局变量k（热扩散率）与内层循环变量k重名，会覆盖热扩散率的初始值，导致贝塞尔函数的参数计算完全错误；
   • 贝塞尔函数应该在每次循环中，根据当前的l、循环变量（建议改名避免冲突）和thau实时计算，而非提前定义。

3. 变量命名冲突
   内层循环使用k作为循环变量，覆盖了全局定义的热扩散率k=1，后续所有涉及热扩散率的计算都会出错，必须将循环变量改为其他名称（如m）。

4. 积分与求和逻辑不符理论

   • 积分函数integer1仅返回1/(sigma**2)，和TPS法中包含贝塞尔函数的积分公式不符，需要对照文献修正积分项；
   • 求和循环的范围range(L,n)和range(1,n)需要与文献中的级数求和范围一致，当前L=1的含义不明确，需确认起始项。

5. 参数定义问题
   λ = 0.169是乙醇的热导率参考值，但如果代码逻辑错误，这个参数无法起到校准作用，需先修正核心逻辑后再验证该值是否合理。

内容的提问来源于stack exchange，提问作者Rosa