瞬时平衡法计算固体火箭发动机内弹道曲线

!计?算?星?型æ装痢?药?的?几?何?尺?寸?

!参?数簓符?号?说μ明é

!-----------------------------------------------------------------------------------!n表括?示?星?角?数簓,num表括?示?将?推?进?剂á沿?肉╝厚?方?向ç分?为a几?等台?分?,m表括?示?选?择?压1力 值μ大洙?小?

!d表括?示?外猘径?,len表括?示?长¤度â,thet表括?示?星?边?夹D角?,epsilon表括?示?角?度â系μ数簓

!r表括?示?过y度â圆2弧?半?径?,r1表括?示?星?角?圆2弧?半?径?,R0表括?示?通 ?用?气?体?常￡数簓

!l表括?示?药?柱ê的?特?征é长¤度â,y0表括?示?初?始?特?征é参?数簓,y1表括?示?燃?尽?特?征é参?数簓

!I0表括?示?总哩?冲?,F表括?示?推?力 ,Poc表括?示?燃?烧?室酣?的?工¤作痢?压1力

!Isp表括?示?比括?冲?,density_p表括?示?密ì度â,k表括?示?比括?热â?比括?rspeed表括?示?燃?速ê，?pn表括?示?压1力 指?数簓

!mpeff表括?示?有瓺效§装痢?药?量?,Cf表括?示?推?力 系μ数簓,Ctz表括?示?特?征é速ê度â,At表括?示?喉æ部?面?积y,a表括?示?燃?速ê系μ数簓

!S表括?示?平?均ê燃?烧?面?积y,e1表括?示?平?均ê肉╝厚?,epsilon1表括?示?减?面?比括?epsilon2表括?示?增?面?比括?

!foresmax表括?示?前©段?最?大洙?相à对?周ì边?长¤,backsmax表括?示?后è段?最?大洙?相à对?周ì边?长¤

!smin表括?示?最?小?相à对?周ì边?长¤,thet1表括?示?周ì边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?

!Ap表括?示?初?始?通 ?气?面?积y,J表括?示?初?始?通 ?气?参?量?,eta表括?示?装痢?填?系μ数簓,Af表括?示?剩骸?药?面?积y,etaf表括?示?剩骸?药?系μ数簓

!de表括?示?每?一?份 肉╝厚?的?长¤度â,Sa表括?示?燃?烧?面?积y数簓组哩?Apa表括?示?通 ?气?面?积y数簓组哩?

!------------------------------------------------------------------------------------program main

implicit none

real(kind=8),parameter :: Pi=3.14

integer :: n,num,i,m

real(kind=8) :: d,len,thet,epsilon,r,r1

real(kind=8) :: l,y0,y1

real(kind=8) :: I0,F,Poc,Pe,R0

real(kind=8) :: Isp,density_p,k,rspeed,Pn

real(kind=8) :: mpeff,Cf,Ctz,At,a

real(kind=8) :: S,e1,epsilon1,epsilon2

real(kind=8) :: foresmax,backsmax,smin,thet1

real(kind=8) :: Ap,J,eta,Af,etaf

real(kind=8) :: error,de

real(kind=8),allocatable :: Sa(:),Apa(:)

!读á入?所ê需â要癮所ê用?参?数簓值μ

open(3,file="design_parameter.dat")

read(3,*) I0 !总哩?冲?

read(3,*) F !推?力

read(3,*) Poc !燃?烧?室酣?压1力

read(3,*) Isp !比括?冲?

read(3,*) density_p !推?进?剂á的?密ì度â

read(3,*) k !比括?热â?比括?

read(3,*) a !燃?速ê系μ数簓

read(3,*) Ctz !读á取?特?征é速ê度âCtz

read(3,*) Pn !压1力 指?数簓

read(3,*) Pe !读á取?喷?管ì出?口ë的?压1力 Pe

read(3,*) R0 !读á取?通 ?用?气?体?常￡数簓R0

read(3,*) d,r !读á取?装痢?药?直ª径?d和æ过y度â圆2弧?半?径?r

read(3,*) epsilon2 !读á取?增?面?比括╡psilon2

read(3,*) r1 !读á取?星?角?圆2弧?半?径?r1

read(3,*) thet1 !试?取?周ì边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?thet1

read(3,*) thet !试?取?初?始?时骸?的?星?边?夹D角?thet

read(3,*) n !读á取?星?角?数簓n（辍?,4,5,6,7,8）?

read(3,*) epsilon !试?取?角?度â系μ数簓epsilon

read(3,*) num !读á取?等台?分?肉╝厚?的?等台?分?数簓num

read(3,*) m !m=1表括?示?最?小?压1力 给?定¨,m=2表括?示?最?大洙?压1力 给?定¨,m=0表括?示?平?均ê压1力 给?定¨

close(3)

allocate(Sa(0:num),Apa(0:num))

!根ê据Y规?定¨的?总哩?冲?计?算?有瓺效§装痢?药?量?mpeff

mpeff=1.02*I0/Isp

!计?算?推?力 系μ数簓和æ喉æ部?面?积y

Cf=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)*&

sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))

At=F/Cf/(Poc*101325.0)

!计?算?平?均ê肉╝厚?e1和æ计?算?特?征é长¤度â和æ平?均ê燃?烧?面?积yS

!若?给?定¨的?是?最?大洙?压1力 则ç需â根ê据Y增?面?比括?计?算?最?小?燃?面?再ê计?算?平?均ê燃?面?值μ

!y1一?般?取?值μ在ë1附?近ì,(0.8-1.2)

if(m==0) then!平?均ê压1力 给?定¨

S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)

e1=mpeff/(density_p*S)

else if(m==1) then!最?小?压1力 给?定¨

S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)

S=(epsilon2+1.0)*S/2.0

e1=mpeff/(density_p*S)