|
|
本帖最后由 xiao. 于 2021-2-4 17:58 编辑
不知道为啥那个变量的值一直是NAN,打断点找原因也找不出来，求助大家！
大噶们帮忙试一试，谢谢啦，再次谢谢大家！
ipl_r 变量在第96行左右！
[Fortran] 纯文本查看 复制代码
!***********************参数设置*************************************************************
implicit real*8(i-n)
real*8 x,a,b,answer
real*8 B_NUM(2966296),AB_CO(2966296) !波数，吸收系数
real*8 fk(10000) !k的概率分布函数
real*8 fk_0 !记录吸收系数的对数小于-22的概率分布函数
real*8 gk(10000) !k的累积概率分布函数
real*8 fki(10000) !kIg的概率分布函数
real*8 gki(10000) !kIg的累积概率分布函数
REAL*8 Ig !g空间上的光谱辐射强度Ig=Ik/f(T,k),Ik=∫Iηδ(k-kη)dη,cf(T,k)=1/Ib(∫Ibηδ(k-kη)dη)
REAL*8 PL_IB !黑体总辐射强度
REAL*8 PL_IBG(2966296) !黑体光谱辐射强度
REAL*8,PARAMETER:: PI=3.141592653589793
REAL*8,PARAMETER:: C1=3.742*10.0**(-8) !普朗克第一常数,(W·cm^4/m^2)
REAL*8,PARAMETER:: C2=1.4388 !普朗克第二常数,(cm·K)
REAL*8,PARAMETER:: E=2.718281828 !自然常数
REAL*8,PARAMETER:: tg=1800.0 !气体温度,K
REAL*8 DELTAK !吸收系数间隔
real*8 maxk !吸收系数最大值
real*8 gaussh !剖分步长
real*8 sum_gk ! 利用累计k分布直接计算k的积分值
real*8 sum_lbl !利用lbl直接计算k的积分值
real*8 sum_gauss !利用高斯求积计算k的积分值
real*8 ipl_r(21) !第i个辐射源项到第j个辐射热流点的辐射热流，从左到右
real*8 ipr_l(21) !第i个辐射源项到第j个辐射热流点的辐射热流，从右到左
real*8 delta_q(20) !辐射源项
real*8 delta_x !控制体的长度,delta_x/划分区间数,cm
real*8 distance !总长度,cm
real*8 temp_sum_l(21),temp_sum_r(21),temp(21)
character(10) line_distance
!**************************************************************
!********外部声明函数******************************************
external expi3 !3阶指数积分函数
!**************************************************************
!*****将数据读取到B_NUM和AB_CO中,并将吸收吸收转为以10为底的对数
OPEN(11,FILE='T_1000_P_1.0_CO2Mol_0.1_H2OMol_0.2_COMol_0.03_N2mol
&_0.67_soot.data')
do I=1,2966296
READ(11,*) B_NUM(I),AB_CO(I)
AB_CO(I)=LOG10(AB_CO(I))
enddo
CLOSE(11)
!*******************计算吸收系数最大值********************************
maxk=0
do i=1,2966295
if(AB_CO(I+1).gt.AB_CO(I).and.AB_CO(I+1).gt.maxk) then
maxk=AB_CO(I+1)
endif
enddo
!*******************计算吸收系数间隔**********************************
DELTAK=(maxk+22.0)/10000.0
!**********黑体总辐射强度,温度为2000K*********************************
PL_IB=(5.67*(10.0**-8.0)*tg**4.0)/PI
!**********黑体光谱辐射强度,单位为(W/(m^2·cm·Sr))*******************
PL_IBG(1:2966296)=0
DO I=1,2966296
PL_IBG(I)=(B_NUM(I)**3.0)*C1/(E**(C2*B_NUM(I)/tg)-1.0)/PI
enddo
!**********fk*dk，记录△log10(k)内的点数******************************
fk(1:10000)=0 !数组初始化为0
fk_0=0 !初始化为0
DO I=1,2966296
if(AB_CO(I).LE.-22.0) then
fk_0=fk_0+0.01*PL_IBG(I)/PL_IB
endif
do J=1,10000
IF(AB_CO(I).GT.(-22.0+(J-1)*DELTAK)
& .AND.AB_CO(I).LT.(-22.0+J*DELTAK) ) THEN
fk(J)=fk(J)+0.01*PL_IBG(I)/PL_IB
ENDIF
enddo
enddo
!**************************************************************************************************
!***********************1维气体辐射传递计算,不考虑壁面辐射*****************************************
!**************************************************************************************************
ipl_r(1:21)=0 !第i个辐射热流点的辐射热流，从左到右
ipr_l(1:21)=0 !第i个辐射热流点的辐射热流，从右到左
delta_q(1:20)=0 !辐射源项
distance=0.1 !总长度,cm
write(line_distance,"(f10.1)") distance
delta_x=distance/20.0 !控制体的长度,delta_x/划分区间数,0.5cm
temp_sum_l(1:21)=0
temp_sum_r(1:21)=0
temp(1:21)=0
!**************************************************************************
! lbl法计算辐射热流项
!**************************************************************************
do i=2,21 !从左到右计算辐射热流点的值
do k=1,2966296
temp_sum_l(1)=10.0**(AB_CO(k))*delta_x*(i-1)
ipl_r(i)=ipl_r(i)+PL_IBG(k)*
& (1-2.0*expi3(temp_sum_l(1)))*0.01
enddo
enddo
do i=20,1,-1 !从右到左计算辐射热流点的值
do k=1,2966296
temp_sum_r(1)=10.0**(AB_CO(k))*(21-i)*delta_x
ipr_l(i)=ipr_l(i)+PL_IBG(k)*
& (1-2.0*expi3(temp_sum_r(1)))*0.01
enddo
enddo
do i=1,21 !计算净辐射热流
temp(i)=pi*(ipl_r(i)-ipr_l(i))
enddo
do i=1,20 !计算辐射源项
delta_q(i)=(temp(i+1)-temp(i))/delta_x/0.01
enddo
open(31,file='插值节点辐射热源'//trim(adjustl(line_distance))//'
&T_1800_P_25.0_H2OMol_0.1_CO2Mol_0.1_N2mol_0.8-kdist.txt')
do i=1,20
write(31,*) delta_q(i)
enddo
close(31)
open(32,file='插值节点净辐射热流'//trim(adjustl(line_distance))//'
&T_1800_P_25.0_H2OMol_0.1_CO2Mol_0.1_N2mol_0.8-kdist.txt')
do i=1,21
write(32,*) temp(i)
enddo
close(32)
END
!**************************************************
!**********3阶指数积分计算函数expi3(x)*************
!**************************************************
function expi3(x)
real*8 x,aux,y,expi3
if(x.lt.1.d0) goto 1
y=1.d0/x
aux=1.d0-y*(((y+3.377358d0)*y+2.052156d0)*y+2.709479d-1)/
&((((y*1.072553d0+5.716943d0)*y+6.945239d0)*y+2.593888d0)*y
&+2.709496d-1)
expi3=aux*y*dexp(-x)
goto 3
1 if(x.le.-3.d0) goto 2
aux=(((((((7.122452d-7*x-1.766345d-6)*x+2.928433d-5)*x
&-2.335379d-4)*x+1.664156d-3)*x-1.041576d-2)*x+5.555682d-2)
&*x-2.500001d-1)*x+9.999999d-1
expi3=-1.d75
if(x.ne.0.d0) expi3=x*aux-dlog(dabs(x))-5.772157d-1
goto 3
2 if(x.gt.-9.d0) then
aux=1.d0-((((5.176245d-2*x+3.061037d0)*x+3.243665d1)*x
&+2.244234d2)*x+2.486697d2)/((((x+3.995161d0)*x+3.893944d1)
&*x+2.263618d1)*x+1.807837d2)
else
y=9.d0/x
aux=1.d0-y*(((y+7.659824d-1)*y-7.271015d-1)*y-1.080693d0)
&/((((y*2.518750d0+1.122927d1)*y+5.921405d0)*y
&-8.666702d0)*y-9.724216d0)
endif
expi3=aux*dexp(-x)/x
3 continue
if(x.eq.0.d0) then
expi3=0.5d0
else
expi3=0.5d0*((1.d0-x)*dexp(-x)+x**2*expi3)
endif
return
end function
|
|
F 币
贡献
权杖
窥视卡
发表于 2021-2-4 17:49:24
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
发表于 2021-2-4 21:27:35
发表于 2021-2-5 08:40:59
楼主