问题描述
我一直在尝试使用 Spyder 上的 Python 复制 O'Dwyer 的论文“具有低维接触的热载流子太阳能电池中的电子和热传输”(Link to O'Dwyer Paper) 中的图 1 和图 2。>
要复制的数字
图 1：w = 1e-5
图 2 = w = 1e-2
方法
要找到吸收器温度 T_H，需要将由于辐射引起的净传入能量流 Qrad 和从热吸收器储层流出的净热流 Qabs 相等. 他们的方程如下：
图 1 和图 2 中的粗线图表示 Wurfel 的解，由以下方程给出：
我在复制图 2 时取得了一些成功，其中 w=1e-2（我的结果如下所示）但没有成功获得图 1，其中 w=1e-5（下面的点和 num_T 指的是绘图点的数量和分别迭代的温度数）。
当 w=1e-2，点数 = 21，num_T = 300 时，我在图 2 中的尝试
我想我目前在尝试使用 w=1e-5 使图 1 正常工作时遇到“exp 中遇到溢出”警告的问题。当我尝试计算 Qabs（请参阅下面“参数”函数中的代码）时，它给出了数量级约为 1e-70 的荒谬值。然而，当我在 WolframAlpha 中运行相同的方程时，我得到了更合理的结果。
例如，当 W = 1e-5、N = 1e12 和电压 = 0 V 时，T_H 值为 ~T_H = 1448K（参见图 1，左上图）。使用 WolframAlpha，我得到 4.54986×10^22对于 Qrad 和 4.83602×10^22 对于 Qabs (WolframAlpha solution for Qrad at w=1e-5,N=1e12,V=0) 和 WolframAlpha solution for Qabs at w=1e-5,V=0))，这是我在 Python 中想要的结果。在下面找到我所有的代码。
所有代码
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider,Button
import matplotlib.ticker as ticker
from scipy.integrate import quad
from scipy.special import expit
import time
from sympy import symbols,Eq,solve
# import warnings
# warnings.filterwarnings("ignore")
t0= time.perf_counter()
directory = r'C:\Users\gyanj\Documents\GADGET BACKUP\University\5th Year\Thesis\Python Simul\ODwyer\Plots'
os.chdir(directory)
c = 3e8 #speed of light,m/s
q = 1.602e-19 # charge of electron,C
h = 6.626e-34/q #Planck's Constant,eVs
k = 8.617e-5 # Boltzmann's Constant,eVK^-1
stefan = 5.67e-8 #Stefan-Boltzmann's Constant,Wm^-2K^-4
T_C = 300 #Cold Reservoir Temperature,K
T_S = 6000 #Sun Temperature,K
Omega = np.pi #Absorption/Emission Solid Angle,sr
A = 1e-4 #Absorber Area,m^2
points = 21 # Number of plotting points
num_T = 300 #Number of temperatures to iterate through
Temperatures = np.linspace(T_C,T_S,num_T) # array of temperatures
E_u = 1 #Average electrochemical potential of system,eV
V = np.linspace(0,1,points) #V applied symetrically across device
max_lim = np.inf# integral upper limit
W = [1e-2] #Transmission function width
N = [1e9,1e10,1e12] #Number of contacts
#Following block used for progress bar (not relevant to calculations)
global total
total = len(W)*len(N)*(points)*len(Temperatures)
progress = 0
counter = 0
full_time = 0
#Object containing all relevant parameters
class param:
def __init__(self,TH,I,P,n,Qrad,Qabs):
self.TH = TH #Hot reservoir/Absorber Temperature,K
self.I = I # Current,A/m^2
self.P = P #Power,W/m^2
self.n = n #Efficiency
self.Qrad = Qrad #net incoming energy flow due to radiation
self.Qabs = Qabs #net heat current flowing out of the hot absorber reservoir
Data = np.empty([len(W),len(N),points],dtype = object) #Contain all param objects
datafile = 'ODwyer.dat'
fout = open(datafile,'w')
fout.write('')
fout.close()
for i in range(len(W)):
for j in range(len(N)):
for x in range(points):
Data[i][j][x] = param(0,0)
# Function Paramaters calculates Qrad,Qabs and I for a given T_H,u_H,u_C,N_contact,w,voltage
def Parameters (T_H,voltage):
eqn1 = lambda E: ((E)**3/(np.exp(E/(k*T_S))-1)-(E)**3/(np.exp(E/(k*T_H))-1))
Qrad = ((2*Omega*A*q)/((h**3)*(c**2)))*quad(eqn1,max_lim)[0]
eqn2 = lambda E:(E-u_H)*(expit(-(E-u_H)/(k*T_H))-expit(-(E-u_C)/(k*T_C)))*(np.exp(-(E-E_u/2)**2/(w)))
Qabs = ((4*N_contact*q)/h)*quad(eqn2,max_lim)[0]
if Qabs < 0:
Qabs = np.inf
error = abs(Qrad-Qabs)
eqn3 = lambda E:(expit(-(E-u_H)/(k*T_H))-expit(-(E-u_C)/(k*T_C)))*(np.exp(-(E-E_u/2)**2/(w)))
I = -((2*N_contact*q)/h)*quad(eqn3,max_lim)[0]/A
fout = open(datafile,'a')
fout.write('%.2e\t%.2e\t%.1f\t%.2f\t%.2e\t%.2e\n'%(w,T_H,voltage,Qabs))
fout.close()
return error,Qabs
#Progress bar for simulation time (not relevant for calculations)
def progressbar(progress):
if (progress >= 0.01):
t1 = time.perf_counter() - t0
full_time = t1*1/progress*100
timeleft = full_time-t1
if timeleft >= 3600:
timelefthrs = int(round(timeleft/3600,0))
timeleftmins = int((timeleft-timelefthrs*3600)%60)
print('\rSimulation Progress: %.2f%%\t Estimated Time Left: %dh %dm '%(progress,timelefthrs,timeleftmins),end='')
elif timeleft >= 60 and timeleft <3600: # in mins
timeleftmins = int(round(timeleft/60,0))
timeleftsecs = int((timeleft-timeleftmins*60)%60)
print('\rSimulation Progress: %.2f%%\t Estimated Time Left: %dm %ds '%(progress,timeleftmins,timeleftsecs),end='')
else:
print('\rSimulation Progress: %.2f%%\t Estimated Time Left: %ds '%(progress,timeleft),end='')
else:
print('\rSimulation Progress: %.2f%%'%(progress),end='')
def Odwyer(index,counter):
for j in range(len(N)):
for i in range(points): #per V
u_H = E_u+V[i]/2 #Hot absorber electrochemical potential,eV
u_C = E_u-V[i]/2 #Cold Reservoir electrochemical potential,eV
error = np.inf #initialise error between Qrad and Qabs as inf
for x in range(len(Temperatures)):
temperature = Temperatures[x]
diff,Qabs= Parameters(Temperatures[x],N[j],W[index],V[i])
if diff <= error: #if difference between Qabs and Qrad is smaller than prevIoUs error,use this Temperature[x]
Data[index][j][i].TH = temperature
Data[index][j][i].Qrad = Qrad
Data[index][j][i].Qabs = Qabs
Data[index][j][i].I = I
Data[index][j][i].P = I*V[i]
Data[index][j][i].n = I*V[i]/(stefan*(T_S**4))
error = abs(diff)
counter += 1
progress = counter/total*100
progressbar(progress)
#Plotting
fig,axs= plt.subplots(2,2,constrained_layout=True)
ax1 = axs[0,0]
ax2 = axs[0,1]
ax3 = axs[1,0]
ax4 = axs[1,1]
for i in range(2):
for j in range(2):
axs[i,j].set_xlim(0,1)
axs[i,j].xaxis.set_major_locator(ticker.MultipleLocator(0.5))
axs[i,j].set_xlabel("Voltage (V)")
ax1.set_ylim(0,T_S)
ax1.set_ylabel("TH (K)")
ax1.yaxis.set_major_locator(ticker.MultipleLocator(2000))
ax2.set_ylim(0,1e8)
ax2.set_ylabel("I (A/m^2)")
ax2.yaxis.set_major_locator(ticker.MultipleLocator(2e7))
ax3.set_ylim(0,1e8)
ax3.set_ylabel("Power (W/m^2)")
ax3.yaxis.set_major_locator(ticker.MultipleLocator(2e7))
ax4.set_ylim(0,1)
ax4.set_ylabel("Efficiency")
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.2))
TH = np.empty([len(N),points])
I = np.empty([len(N),points])
P = np.empty([len(N),points])
n = np.empty([len(N),points])
for j in range(len(N)):
for x in range(points):
TH[j][x] = Data[index][j][x].TH
I[j][x] = Data[index][j][x].I
P[j][x] = Data[index][j][x].P
n[j][x] = Data[index][j][x].n
#Wurfel's Solution
TH_W = []
I_W = []
P_W = []
n_W = []
for x in range(points):
if V[x] == E_u:
TH_wurfel = 1e20
else:
TH_wurfel = T_C/(1-V[x]/E_u)
TH_W.append(TH_wurfel)
Iwurfel = (stefan)/(E_u)*(T_S**4-TH_wurfel**4)
Pwurfel = stefan*(T_S**4-TH_wurfel**4)*(1-T_C/TH_wurfel)
nwurfel = (T_S**4-TH_wurfel**4)/(T_S**4)*(1-T_C/TH_wurfel)
I_W.append(Iwurfel)
P_W.append(Pwurfel)
n_W.append(nwurfel)
linestyles = ['--','-','-.']
for j in range(len(N)):
for x in range(points):
if TH[j][x] == T_S:
TH[j][x] = 1e8
for i in range(len(N)):
ax1.plot(V,TH[i],label='N = %.0e'%N[i],color = 'black',linestyle = linestyles[i],linewidth = 1)
ax2.plot(V,I[i],linewidth = 1)
ax3.plot(V,P[i],linewidth = 1)
ax4.plot(V,n[i],linewidth = 1)
ax1.plot(V,TH_W,label='Wurfel',linewidth = 3)
ax2.plot(V,I_W,linewidth = 3)
ax3.plot(V,P_W,linewidth = 3)
ax4.plot(V,n_W,linewidth = 3)
fig.suptitle('w = %.0e eV' % W[index])
ax1.legend(loc='upper right',fontsize = 8)
ax2.legend(loc='upper right',fontsize = 8)
ax3.legend(loc='upper right',fontsize = 8)
ax4.legend(loc='upper right',fontsize = 8)
#Saving figure
fig.savefig('w = %.0e eV,pp = %d,num_T = %d.jpg' %(W[index],points,num_T),dpi=800)
return counter
for x in range(len(W)):
counter = Odwyer(x,counter)
# Printing out object values
for x in range(len(W)):
for j in range(len(N)):
print('Parameters for W = %0.e,N = %.0e'%(W[x],N[j]))
for i in range(points):
print('w = %.0e\tV = %.2f\tTH = %.0f\tQrad = %.2e\tQabs = %.2e\tI = %.2e'%(W[x],V[i],Data[x][j][i].TH,Data[x][j][i].Qrad,Data[x][j][i].Qabs,Data[x][j][i].I))
print('\nComplete!')
我的尝试
我尝试将积分的上限从 inf 更改为较低的值，尽管它删除了值 ~Overflow Post 等溢出的解决方案，但无济于事。这是我的第一篇文章，所以如果您需要我提供更多信息来解决我的问题，请随时告诉我！
解决方法
这是一个快速而肮脏的 stdlib（没有 numpy）脚本，它与 WolframAlpha 的答案很接近：
from math import exp,pi
C1 = 8.617e-5 * 6000
C2 = 8.617e-5 * 1448
def f(x):
denom1 = exp(x / C1)
denom2 = exp(x / C2)
# did some algebra
difference = (denom2 - denom1) / (denom1 - 1) / (denom2 - 1)
return x ** 3 * difference
bins = 10_000
endpoint = 10
total = 0.0
for i in range(1,bins+1):
x = i * endpoint / bins
total += f(x)
# account for widths
total *= (endpoint / bins)
scaled = float(total) * 2 * pi * 1e-4 / (4.14e-15)**3 / (3e8)**2
print(scaled)
# 4.549838698077388e+22
问题的一部分（我猜，不确定）是 1/(a-1) - 1/(b-1) 对于接近 1 的 a 和 b 会非常不精确，所以你可以做一些代数来尝试解决这个问题，并使它(b-a)/(a-1)/(b-1)。