本程序从固体火箭发动机假定化学式计算到燃烧室热力计算,再到喷
管热力计算,最后到星形药柱几何参量计算和内弹道计算。
!从推进剂假定化学式开始到喷气速度的计算
!计?算?燃?烧?室酣?和í喷?管ü过y程ì中D一?些?热è?力 参?量?
!固ì体?火e箭y发ぁ?动ˉ机ú燃?烧?室酣?喷?管ü的?热è?力 计?算?程ì序ò
!计?算?推?进?剂á装痢?药?的?几?何?参?数簓值μ
!计?算?相à应畖装痢?药?的?内ú弹獭?道台?曲ú线?
!全?部?使?用?国ú际ê单蹋?位?制?
!-----------------------------------------------------------------------------
!L表示组分凝相数,M表示推进剂化学元素,N表示总组分数,KK表示组元数
!Roc表括?示?燃?烧?室酣?的?燃?烧?产ú物?常£数簓(J/Kg*K),Soc表括?示?燃?烧?室酣?的?总哩?
熵?值μ
!MM表括?示?喷?管ü出?口ú不?考?虑?的?次?要癮组哩?分?数簓,std=1为a选?择?平?衡a流ⅰ?否?则
ò为a冻3结á流ⅰ?
!Pe为a喷?管ü出?口ú压1力 ,Rme为a出?口ú气?体?常£数簓(J/Kg*K),k为a等台?熵?指?数簓,Ue为
a出?口ú气?流ⅰ?速ù度è
!Ime为a出?口ú总哩?焓 ?值μ,Te为a出?口ú温?度è,Xe表括?示?喷?管ü出?口ú气?体?的?组哩?分?
摩|尔?数簓
!Ip表括?示?推?进?剂á在ú初?温?下?的?总哩?焓 ?Toc表括?示?燃?烧?室酣?实害?际ê的?温?度è
值μ
!Im1表括?示?燃?烧?温?度è为aT1时骸?组哩?分?的?总哩?焓 ?值μ,Im2表括?示?燃?烧?温?度è为
aT2时骸?组哩?分?的?总哩?焓 ?值μ
!Ctz表括?示?燃?烧?室酣?理え?论?特?征÷速ù度è,Cp表括?示?定¨压1比括?热è?Cv表括?示?定¨
容╕比括?热è?Moc表括?示?燃?烧?室酣?燃?气?分?子哩?量?
!Moc1表括?示?燃?烧?室酣?内ú凝y相à组哩?分?的?分?子哩?量?,percent表括?示?凝y相à组哩?分?
的?质ê量?百悒?分?数簓
!A1计?算?燃?烧?室酣?焓 ?值μ和í熵?值μ代洙?数簓式?中D的?系μ数簓矩?阵ó
!EX表括?示?组哩?元a的?百悒?分?数簓,EC表括?示?组哩?元a一?般?化ˉ学§式?,EA表括?示?元a
素?原-子哩?数簓
!----------------------------------------------------------------------------
program main
implicit none
integer :: i,j,L,M,N
integer :: KK
real(kind=8) :: Roc,Soc
integer :: MM,std
real(kind=8) :: Pe,Rme,k,Ue,Ime,Te
real(kind=8),allocatable :: Xe(:)
character(len=8),allocatable :: string1(:),string2(:),string(:)
real(kind=8) :: R0,Toc,Poc,Ip,error
real(kind=8) :: T1,Im1,T2,Im2,Xg
real(kind=8) :: Ctz,Cp,Cv,Moc,percent
real(kind=8),allocatable :: Moc1(:)
real(kind=8),allocatable :: NK(:),A(:,:),C(:),X(:),X1(:),X2(:),Cp1(:)
real(kind=8),allocatable :: H1(:),H0(:),S1(:),H2(:),S2(:),H(:),S(:)
real(kind=8),allocatable :: A1(:,:)
real(kind=8),allocatable :: EX(:),EC(:,:),EA(:)
!从洙?文?件t夹D中D读á入?参?数簓
open(3,file="parameter.dat")
read(3,*) L,M,N,KK !读á入?凝y相àL,推?进?剂á组哩?成é元a素?M,组哩?分?数簓N,KK表括?示?组
哩?元a数簓
!读á入?喷?管ü参?数簓MM表括?示?不?考?虑?的?次?要癮组哩?分?数簓,std=1为a平?衡a流ⅰ?否?
则ò为a冻3结á流ⅰ?k为a试?算?的?绝?热è?指?数簓,Pe为a出?口ú压1力
read(3,*) MM,std,k,Pe
read(3,*) R0,Ip,Poc
!读á入?通 ?用?气?体?常£数簓R0,推?进?剂á总哩?焓 ?值μIp,燃?烧?
室酣?压1力 Poc
read(3,*) T1,T2 !输?入?第台?一?次?计?算?燃?烧?室酣?组哩?分?的?两?个?试?算?温?度è
read(3,*) error !读á入?收?敛?误ó差?error
close(3)
allocate(string1(N)) !表括?示?燃?烧?室酣?燃?烧?产ú物?的?化ˉ学§式?
allocate(string2(M)) !表括?示?推?进?剂á元a素?的?化ˉ学§符?号?
allocate(string(N-MM)) !表括?示?喷?管ü燃?烧?产ú物?的?化ˉ学§式?
allocate(Cp1(N)) !表括?示?燃?烧?产ú物?的?定¨压1比括?热è?
allocate(Moc1(L)) !表括?示?燃?烧?室酣?组哩?分?中D凝y相à的?摩|尔?分?子哩?量?
allocate(C(N)) !表括?示?迭台?代洙?前°一?次?的?组哩?分?摩|尔?数簓值μ
allocate(X(N)) !表括?示?温?度èTf时骸?这a次?迭台?代洙?计?算?的?组哩?分?摩|尔?数簓值μ
allocate(X1(N)) !表括?示?温?度èT1时骸?这a次?迭台?代洙?计?算?的?组哩?分?摩|尔?数簓值μ
allocate(X2(N)) !表括?示?温?度èT2时骸?这a次?迭台?代洙?计?算?的?组哩?分?摩|尔?数簓值μ
allocate(Xe(N-MM)) !表括?示?喷?管ü出?口ú组哩?分?的?摩|尔?数簓值μ
allocate(NK(M)) !表括?示?推?进?剂á假ù定¨化ˉ学§式?中D元a素?的?摩|尔?原-子哩?数簓
allocate(A(M,N)) !表括?示?某3种?组哩?分?中D某3种?元a素?的?原-子哩?数簓
allocate(H1(N)) !表括?示?T1温?度è燃?烧?组哩?分?的?总哩?焓 ?值μ
allocate(H0(N)) !表括?示?298.15K下?的?标括?准?摩|尔?生Θ?成é焓 ?
allocate(S1(N)) !表括?示?T1温?度è组哩?分?在ú一?个?物?理え?大洙?气?压1下?的?熵?值μ
allocate(H2(N)) !表括?示?T2温?度è燃?烧?组哩?分?的?总哩?焓 ?值μ
allocate(S2(N)) !表括?示?T2温?度è组哩?分?在ú一?个?物?理え?大洙?气?压1下?的?熵?值μ
allocate(H(N)) !表括?示?Tf温?度è燃?烧?组哩?分?的?总哩?焓 ?值μ
allocate(S(N)) !表括?示?Tf温?度è组哩?分?在ú一?个?物?理え?大洙?气?压1下?的?熵?值μ
allocate(A1(N,8)) !表括?示?燃?烧?室酣?计?算?焓 ?值μ和í熵?值μ的?系μ数簓矩?阵ó
allocate(EX(KK)) !表括?示?组哩?元a百悒?分?数簓
allocate(EC(M,KK)) !表括?示?组哩?元a一?般?化ˉ学§式?
allocate(EA(M)) !表括?示?元a素?的?原-子哩?量?
!从洙?文?件t夹D中D读á入?燃?烧?组哩?分?的?摩|尔?原-子哩?数簓A
open(7,file="A.dat")
do i=1,M,1
read(7,*) A(i,1:N)
end do
close(7)
!从洙?文?件t中D读á取?试?算?值μC和í化ˉ学§式?
open(8,file="stringC.dat")
read(8,*) string1(1:N)
read(8,*) C(1:N)
close(8)
!读á取?燃?烧?室酣?燃?烧?温?度èT时骸?计?算?焓 ?值μ和í熵?值μ的?系μ数簓矩?阵ó
open(9,file="A_oc.dat")
do i=1,N,1
read(9,*) A1(i,1:8)
end do
close(9)
!读á取?298.15K下?的?标括?准?摩|尔?生Θ?成é焓 ?
open(11,file="H0.dat")
read(11,*) H0(1:N)
close(11)
!读á取?组哩?元a百悒?分?数簓EX和í元a素?原-子哩?量?EA的?值μ
open(12,file="EAEX.dat")
read(12,*) string2(1:M)
read(12,*) EA(1:M)
read(12,*) EX(1:KK)
close(12)
!读á取?组哩?元a的?一?般?化ˉ学§式?EC
open(13,file="EC.dat")
do i=1,M,1
read(13,*) EC(i,1:KK)
end do
close(13)
write(*,*) "请?耐í心?等台?待鋣程ì序ò执′行D!?"
!调獭?用?假ù定¨化ˉ学§式?求ó解a程ì序ò计?算?NK的?值μ
call formula(KK,M,EX,EC,EA,NK)
open(10,file="exportNK.dat")
do i=1,M,1
write(10,"(A2,2x,f9.5)") string2(i),NK(i)
end do
close(10)
!循-环·计?算?燃?烧?室酣?推?进?剂á燃?烧?时骸?的?组哩?分?和í温?度è,?直±到?满ú足?要癮
求ó为a止1
do while(.true.)
!计?算?温?度è为aT1时骸?的?总哩?焓 ?值μ(J/mol)和í一?个?物?理え?大洙?气?压1熵?值μ
(J/mol*K)
do i=1,N,1
H1(i)=1000.0*(A1(i,1)*(T1/1000.0)+A1(i,2)*(T1/1000.0)**2/2.0+A1(i,3)*(T1/1000.0)**3/3+&
A1(i,4)*(T1/1000.0)**4/4-A1(i,5)*1000.0/T1+A1(i,6)-A1(i,8)+H0(i))
S1(i)=A1(i,1)*log(T1/1000.0)+A1(i,2)*T1/1000.0+A1(i,3)*(T1/1000.0)**2/2.0+&
A1(i,4)*(T1/1000.0)**3/3-A1(i,5)/(2*(T1/1000.0)**2)+A1(i,7)
end do
!计?算?温?度è为aT2时骸?的?总哩?焓 ?值μ(J/mol)和í一?个?物?理え?大洙?气?压1熵?值μ
(J/mol*K)
do i=1,N,1
H2(i)=1000.0*(A1(i,1)*(T2/1000.0)+A1(i,2)*(T2/1000.0)**2/2.0+A1(i,3)*(T2/1000.0)**3/3+&
A1(i,4)*(T2/1000.0)**4/4-A1(i,5)*1000.0/T2+A1(i,6)-A1(i,8)+H0(i))
S2(i)=A1(i,1)*log(T2/1000.0)+A1(i,2)*T2/1000.0+A1(i,3)*(T2/1000.0)**2/2.0+&
A1(i,4)*(T2/1000.0)**3/3-A1(i,5)/(2*(T2/1000.0)**2)+A1(i,7)
end do
!调獭?用?最?小?自?由 ?能ü法ぁ?子哩?程ì序ò求ó解aT1温?度è时骸?组哩?分?的?摩|尔?数簓
call mfem(error,R0,T1,Poc,L,M,N,A,C,NK,H1,S1,X1)
!计?算?温?度è为aT1时骸?组哩?分?的?总哩?焓 ?值μIm1
Im1=0.0
do i=1,N,1
Im1=Im1+H1(i)*X1(i)
end do
!调獭?用?最?小?自?由 ?能ü法ぁ?子哩?程ì序ò求ó解aT2温?度è时骸?组哩?分?的?摩|尔?数簓
call mfem(error,R0,T2,Poc,L,M,N,A,C,NK,H2,S2,X2)
!计?算?温?度è为aT1时骸?组哩?分?的?总哩?焓 ?值μIm1
Im2=0.0
do i=1,N,1
Im2=Im2+H2(i)*X2(i)
end do
!判D断?是?否?满ú足?条?件t利?用?牛£顿ù内ú插?值μ法ぁ?计?算?温?度è
if(Ip>Im1.and.Ip
A1(i,4)*(Toc/1000.0)**3/3-A1(i,5)/(2*(Toc/1000.0)**2)+A1(i,7)
end do
!调獭?用?最?小?自?由 ?能ü法ぁ?子哩?程ì序ò求ó解aTf温?度è时骸?组哩?分?的?摩|尔?
数簓
call mfem(error,R0,Toc,Poc,L,M,N,A,C,NK,H,S,X)
exit
else if(Ip==Im1) then
Toc=T1
X=X1
exit
else if(Ip==Im2) then
Toc=T2
X=X2
exit
else if(Ip
!-------------------------------------------------
!计?算?星?型í装痢?药?的?几?何?尺?寸?
!-------------------------------------------------
call xin_ge(Poc,Pe,R0,Ue,Ctz,k,Roc,Toc)
!-------------------------------------------------
!---------------------------------------
!计?算?燃?烧?产ú物?的?热è?力 参?数簓
!---------------------------------------
!计?算?燃?烧?室酣?凝y相à组哩?分?的?质ê量?百悒?分?数簓ε?值μ
Moc1=0.0
do i=1,L,1
do j=1,M,1
Moc1(i)=Moc1(i)+A(j,i)*EA(j)
end do
end do
percent=0.0
do i=1,L,1
percent=percent+Moc1(i)*X(i)/1000.0
end do
!计?算?燃?烧?产ú物?气?相à的?平?均ù摩|尔?质ê量?(Kg/mol)
Moc=(1-percent)/Xg
!各÷组哩?分?的?定¨压1比括?热è?值μCp
do i=1,N,1
Cp1(i)=A1(i,1)+A1(i,2)*Toc/1000.0+A1(i,3)*(Toc/1000.0)**2+&
A1(i,4)*(Toc/1000.0)**3+A1(i,5)/(Toc/1000.0)**2
end do
!计?算?燃?烧?产ú物?的?平?均ù定¨压1比括?热è–p和í定¨容╕比括?热è–v
Cp=0.0
do i=1,N,1
Cp=Cp+X(i)*Cp1(i)
end do
Cv=Cp-R0/Moc
!输?出?燃?烧?室酣?的?组哩?分?、¢气?体?常£数簓、¢总哩?熵?和í燃?烧?温?度è的?值μ
open(20,file="exportX_oc.dat")
do i=1,N,1
write(20,"(A8,2x,f8.5,A6)") string1(i),X(i),"mol/Kg"
end do
write(20,"(A4,2x,f9.5,A6)") "Roc=",Roc,"J/Kg*K"
write(20,"(A4,2x,f10.5,A6)") "Soc=",Soc,"J/Kg*K"
write(20,"(A4,2x,f10.5,A1)") "Toc=",Toc,"K"
write(20,"(A3,2x,f10.5,A3)") "C*=",Ctz,"m/s"
write(20,"(A3,2x,f10.5,A6)") "Cp=",Cp,"J/Kg*K"
write(20,"(A3,2x,f10.5,A6)") "Cv=",Cv,"J/Kg*K"
write(20,"(A4,2x,f8.5,A6)") "Moc=",Moc,"Kg/mol"
close(20)
!输?出?喷?管ü参?数簓
open(30,file="export_e.dat")
do i=1,N-MM,1
write(30,"(A8,2x,f8.5,A6)") string(i),Xe(i),"mol/Kg"
end do
write(30,"(A3,2x,f10.5,A3)") "Ue=",Ue,"m/s"
write(30,"(A4,2x,f9.5,A6)") "Rme=",Rme,"J/Kg*K"
write(30,"(A2,2x,f7.5)") "k=",k
close(30)
stop
end
!该?程ì序ò用?于 ?计?算?推?进?剂á的?假ù定¨化ˉ学§式?
subroutine formula(KK,M,EX,EC,EA,NK)
implicit none
integer :: i,j,KK,M
real(kind=8),intent(in) :: EX(KK),EA(M)
real(kind=8) :: MM(KK),EC(M,KK)
real(kind=8),intent(out) :: NK(M)
!计?算?组哩?元a的?摩|尔?原-子哩?量?
MM=0.0
do i=1,KK,1
do j=1,M,1
MM(i)=MM(i)+EC(j,i)*EA(j)
end do
end do
!计?算?组哩?元a的?假ù定¨化ˉ学§式?
do i=1,KK,1
do j=1,M,1
EC(j,i)=EC(j,i)*1000.0/MM(i)
end do
end do
!计?算?推?进?剂á的?假ù定¨化ˉ学§式?
NK=0.0
do i=1,M,1
do j=1,KK,1
NK(i)=NK(i)+EC(i,j)*EX(j)
end do
end do
return
end subroutine
!高?斯1消?元a法ぁ?
subroutine gauss(n,A,B)
implicit none
integer :: i,j,k,n,l=1,js(n),is
real :: t,d
real(kind=8) :: A(n,n),B(n)
do k=1,n-1,1
d=0.0
do i=k,n,1
do j=k,n,1
if(abs(A(i,j))>d) then
d=abs(A(i,j))
js(k)=j
is=i
end if
end do
end do
if(d+1.0==1.0) then
l=0
else
if(js(k)/=k) then
do i=1,n
t=A(i,k)
A(i,k)=A(i,js(k))
A(i,js(k))=t
end do
end if
if(is/=k) then
do j=k,n
t=A(k,j)
A(k,j)=A(is,j)
A(is,j)=t
end do
t=B(k)
B(k)=B(is)
B(is)=t
end if
end if
if(l==0) then
write(*,"(A6)") "fail"