我一直在尝试复制O'Dwyer论文中的图1和图2；具有低尺寸接触的热载体太阳能电池中的电子和热传输"；(链接到O’Dwyer论文(，在Spyder上使用Python。
要复制的数字
图1:w=1e-5
图1
图2=w=1e-2
图2
方法
为了求出吸收器温度T_H，需要将由于辐射引起的净入射能量流Qrad与流出热吸收器储器的<em]净热流Qabs相等>
Qrad和Qabs 方程
图1和图2中的粗线图指的是由以下方程给出的Wurfel解：
Wurfel的解决方案
我成功地复制了图2，其中w=1e-2(我的结果如下所示(，但没有成功地获得图1，其中w＝1e-5(下面的点和num_T分别指绘制点的数量和要迭代的温度的数量(。
当w=1e-2，points=21，num_T=300时，我在图2中的尝试
我对图2 的尝试
我想我目前在"；在exp"中遇到溢出；警告试图使w=1e-5的图1工作。当我试图计算Qabs(请参阅下面"参数"函数中的代码(时，它给出了数量级为~1e-70的荒谬值。然而，当我在WolframAlpha中运行相同的方程时，我得到了一个更合理的结果。
例如，当W=1e-5，N=1e12，电压=0V时，T_H值为~ T_H=1448K(参见图1左上图(。使用WolframAlpha，我得到了Qrad的4.54986×10^2和Qabs的4.83602×10^2(Qrad在w=1e-5，N=1e12，V=0时的WolframAlpha解(，以及Qabs在w=1e-5，N=1e12，V=1(时的WolframAlpha解，这是我在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:UsersgyanjDocumentsGADGET BACKUPUniversity5th YearThesisPython SimulODwyerPlots'
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,0,0,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, u_H, u_C, N_contact, w, 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,0,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,0,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,0,max_lim)[0]/A
fout = open(datafile,'a')
fout.write('%.2et%.2et%.1ft%.2ft%.2et%.2en'%(w,N_contact,T_H,voltage,Qrad,Qabs))
fout.close()
return error, I, Qrad, 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, I, Qrad, Qabs= Parameters(Temperatures[x], u_H, u_C, 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], label='N = %.0e'%N[i], color = 'black', linestyle = linestyles[i], linewidth = 1)
ax3.plot(V,P[i], label='N = %.0e'%N[i], color = 'black', linestyle = linestyles[i], linewidth = 1)
ax4.plot(V,n[i], label='N = %.0e'%N[i], color = 'black', linestyle = linestyles[i], linewidth = 1)
ax1.plot(V,TH_W, color = 'black', label='Wurfel', linewidth = 3)
ax2.plot(V,I_W, color = 'black', label='Wurfel', linewidth = 3)
ax3.plot(V,P_W, color = 'black', label='Wurfel', linewidth = 3)
ax4.plot(V,n_W, color = 'black', label='Wurfel', 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 = %.0etV = %.2ftTH = %.0ftQrad = %.2etQabs = %.2etI = %.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更改为较低的值，尽管它删除了值的溢出警告~<15，则Qabs=0.00e00。我还尝试在"quad"函数中更改"limit"one_answers"epsabs"的参数，但也无法实现。更改变量"points"one_answers"num_T"也没有提高我的值的准确性。我也阅读并尝试过相关帖子中关于溢出的解决方案，如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
部分问题(我猜，不确定(是，对于a和b接近1的情况，1/(a-1) - 1/(b-1)将非常不精确，所以你可以做一些代数来尝试解决这个问题，并使其成为(b-a)/(a-1)/(b-1)。