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

!计?算?星?型æ装痢?药?的?几?何?尺?寸?!参?数簓符?号?说μ明é!-----------------------------------------------------------------------------------!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 mainimplicit nonereal(kind=8),parameter :: Pi=3.14integer :: n,num,i,mreal(kind=8) :: d,len,thet,epsilon,r,r1real(kind=8) :: l,y0,y1real(kind=8) :: I0,F,Poc,Pe,R0real(kind=8) :: Isp,density_p,k,rspeed,Pnreal(kind=8) :: mpeff,Cf,Ctz,At,areal(kind=8) :: S,e1,epsilon1,epsilon2real(kind=8) :: foresmax,backsmax,smin,thet1real(kind=8) :: Ap,J,eta,Af,etafreal(kind=8) :: error,dereal(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 !读á取?特?征é速ê度âCtzread(3,*) Pn !压1力 指?数簓read(3,*) Pe !读á取?喷?管ì出?口ë的?压1力 Peread(3,*) R0 !读á取?通 ?用?气?体?常￡数簓R0read(3,*) d,r !读á取?装痢?药?直ª径?d和æ过y度â圆2弧?半?径?rread(3,*) epsilon2 !读á取?增?面?比括╡psilon2read(3,*) r1 !读á取?星?角?圆2弧?半?径?r1read(3,*) thet1 !试?取?周ì边?长¤取?得?最?小?值μ时骸?的?星?边?夹D角?thet1read(3,*) thet !试?取?初?始?时骸?的?星?边?夹D角?thetread(3,*) n !读á取?星?角?数簓n（辍?,4,5,6,7,8）?read(3,*) epsilon !试?取?角?度â系μ数簓epsilonread(3,*) num !读á取?等台?分?肉╝厚?的?等台?分?数簓numread(3,*) m !m=1表括?示?最?小?压1力 给?定¨,m=2表括?示?最?大洙?压1力 给?定¨,m=0表括?示?平?均ê压1力 给?定¨close(3)allocate(Sa(0:num),Apa(0:num))!根ê据Y规?定¨的?总哩?冲?计?算?有瓺效§装痢?药?量?mpeffmpeff=1.02*I0/Isp!计?算?推?力 系μ数簓和æ喉æ部?面?积yCf=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.0e1=mpeff/(density_p*S)else!最?大洙?压1力 给?定¨S=At*(Poc*101325.0)**(1-Pn)/(Ctz*density_p*a)S=(1.0/epsilon2+1.0)*S/2.0e1=mpeff/(density_p*S)end ifl=d/2-e1-ry0=r/ly1=(e1+r)/l!根ê据Yepsilon2选?择?角?度â系μ数簓epsilon!计?算?取?得?最?小?周ì长¤时骸?的?thet1do while(.true.)error=thet1/2+cotan(thet1/2)-Pi/n-Pi/2if(error>=0.0001) thenthet1=thet1+0.00001else if(error<=-0.0001) thenthet1=thet1-0.00001elseexitend ifend do!计?算?角?度â系μ数簓epsilon的?准?确í?值μdo while(.true.)backsmax=2*n*((1-epsilon)*Pi/n+y1*(Pi/n+asin(sin(epsilon*Pi/n)/y1))) smin=2*n*(sin(epsilon*Pi/n)/sin(thet1/2)+(1-epsilon)*Pi/n) if(backSmax/Smin>epsilon2) thenepsilon=epsilon+0.001elseexitend ifend do!根ê据Y减?面?比括╡psilon1，?计?算?thetepsilon1=1/epsilon2foresmax=smin/epsilon1do while(.true.)error=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cos(thet/2)/sin(thet/2))/l)-foresmax if(error>=0.0001) thenthet=thet+0.00001else if(error<=-0.0001) thenthet=thet-0.00001elseexitend ifend do!计?算?初?始?通 ?气?面?积y，?通 ?气?参?量?，?装痢?填?系μ数簓,药?柱ê长¤度âAp=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*r*(sin(epsilon*Pi/n)/sin(thet/2)+&(1-epsilon)*Pi/n)/l+n*r**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+&n*r1**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2J=At/Apeta=4*(Pi*d**2/4-Ap)/(Pi*d**2)Af=(epsilon*Pi*(1+y1)**2-n*(sin(epsilon*Pi/n)*(sqrt(y1**2-sin(epsilon*Pi/n)**2)+& cos(epsilon*Pi/n)))-n*y1**2*(epsilon*Pi/n+asin(sin(epsilon*Pi/n)/y1)))*l**2etaf=4*Af/(pi*d**2)len=mpeff/density_p/(Pi*d**2/4-Ap-Af)!计?算?推?进?剂á的?燃?面?变?化ˉ规?律ã并￠输?出?结á果?de=e1/numopen(10,file="export_burnS.dat")open(20,file="export_Ap.dat")do i=0,num,1if((i*de)<=r1) thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(r1+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l-(r1-i*de)*Pi/n/l)*l*len Apa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2+n*(r1-i*de)**2*(thet/2+cotan(thet/2 )-Pi/2)/l**2)*l**2else if((i*de)>r1.and.(i*de)<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r)) thenSa(i)=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(i*de+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa(i)=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-sin(epsilon*Pi/n)* cotan(thet/2)))+&2*n*(r+i*de)*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n)/l+&n*(r+i*de)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2else if((i*de)>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.(i*de)<=e1) thenSa(i)=2*n*((1-epsilon)*Pi/n+(r+i*de)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(i*de+r)))/l)*l *lenApa(i)=n*((1-epsilon)*Pi*(1+(r+i*de)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+i*de)**2/l* *2-&sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+i*de)**2*(epsilon*Pi/n+&asin(l*sin(epsilon*Pi/n)/(r+i*de)))/l**2)*l**2end ifwrite(10,"(f10.5,2X,f15.5)") i*de,Sa(i)write(20,"(f10.5,2X,f15.5)") i*de,Apa(i)end doclose(10)close(20)!调獭?用?子哩?程å序ç计?算?装痢?药?内ë弹獭?道台?曲ë线?call internal_ballistics0(num,de,Sa,Apa,len,a,pn,Ctz,At,k,density_p,Isp,Pe)!输?出?星?型æ装痢?药?的?几?何?参?数簓open(30,file="export_star_geometry.dat")write(30,"(A6,f11.5,A2)") "mpeff=",mpeff,"Kg"write(30,"(A2,f8.5,A2)") "D=",d,"m"write(30,"(A2,f8.5,A2)") "L=",len,"m"write(30,"(A2,I2)") "n=",nwrite(30,"(A3,f8.5)") "θâ=",thetwrite(30,"(A3,f7.5)") "ε?=",epsilonwrite(30,"(A2,f7.5,A2)") "r=",r,"m"write(30,"(A3,f7.5,A2)") "r1=",r1,"m"write(30,"(A3,f8.5,A2)") "e1=",e1,"m"write(30,"(A3,f10.5,A2)") "Ap=",Ap,"㎡O"write(30,"(A2,f10.5)") "J=",Jwrite(30,"(A3,f7.5)") "η?=",etawrite(30,"(A3,f10.5,A2)") "Af=",Af,"㎡O"write(30,"(A4,f10.5)") "η?f=",etafwrite(30,"(A2,f10.5)") "l=",lwrite(30,"(A3,f10.5)") "y0=",y0write(30,"(A3,f10.5)") "y1=",y1close(30)stopend!该?子哩?程å序ç用?于 ?计?算?零?维?变?截?面?燃?烧?装痢?药?的?内ë弹獭?道台?!该?子哩?程å序ç利?用?瞬2时骸?平?衡a法ぁ?计?算?内ë弹獭?道台?!同?时骸?可ã以?计?算?侵?蚀骸?燃?烧?的?影 ?响å!--------------------------------------------------------------------------!Pc_eq表括?示?工¤作痢?段?的?平?衡a压1力 值μ,dPc_eq表括?示?压1力 修T正y值μ,rspeed表括?示?燃?速ê!time表括?示?工¤作痢?段?的?时骸?间?,Toc表括?示?燃?气?温?度â,Roc表括?示?燃?气?气?体?常￡数簓,Pig表括?示?点?火e压1力!P1_eq表括?示?起e始?段?的?压1力 值μ,P2_eq表括?示?拖?尾2段?的?压1力 值μ,Cf表括?示?推?力 系μ数簓!time1表括?示?起e始?段?的?时骸?间?,time2表括?示?拖?尾2段?的?时骸?间?,Isp表括?示?理え?论?比括?冲?值μ!Feq表括?示?工¤作痢?段?的?推?力 值μ,F1表括?示?起e始?段?的?推?力 值μ,F2表括?示?拖?尾2段?的?推?力 值μ!Poc表括?示?燃?烧?室酣?压1力 值μ,Pe表括?示?喷?管ì出?口ë压1力!--------------------------------------------------------------------------subroutine internal_ballistics0(num,de,Sa,Apa,len,a,pn,Ctz,At,k,density_p,Isp,Pe) implicit noneinteger,parameter :: n=10integer :: num,i,keyreal(kind=8),intent(in) :: de,len,a,pn,Ctz,At,k,density_preal(kind=8),intent(in) :: Sa(0:num),Apa(0:num),Isp,Pereal(kind=8) :: Pc_eq(0:num),dPc_eq(0:num),rspeed(0:num)real(kind=8) :: time(0:num),Toc,Roc,Feq(0:num),Cfreal(kind=8),allocatable :: P1_eq(:),P2_eq(:),F1(:)real(kind=8),allocatable :: time1(:),time2(:),F2(:)real(kind=8) :: ga,Pigreal(kind=8),external :: funopen(3,file="inter_parameter.dat")read(3,*) key !读á入?key的?值μread(3,*) Pig !读á入?点?火e压1强?值μread(3,*) Toc !读á入?绝?热â?燃?烧?温?度â值μread(3,*) Roc !读á入?燃?气?的?气?体?常￡数簓值μallocate(P1_eq(n),F1(n),time1(n))allocate(P2_eq(0:n-1),F2(0:n-1),time2(0:n-1))ga=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)Pig=Pig*101325.0!-------------------------------------------------!计?算?工¤作痢?段?的?压1力 和æ时骸?间?值μ!-------------------------------------------------!计?算?各é燃?层?时骸?燃?烧?室酣?的?平?衡a压1力 Pc_eqdo i=0,num,1Pc_eq(i)=(density_p*Ctz*fun(Sa(i),Apa(i),key)*&a*Sa(i)/At)**(1/(1-pn))end do!计?算?各é燃?层?时骸?平?衡a压1力 的?修T正y值μdPc_eqdPc_eq=0.0do i=1,num,1dPc_eq(i)=Pc_eq(i)*(Sa(i)-Sa(i-1))/(Sa(i)*(1-pn))end do!计?算?各é燃?层?的?实害?际ä平?衡a压1力 值μdo i=0,num,1Pc_eq(i)=Pc_eq(i)-Apa(i)*len*fun(Sa(i),Apa(i),key)*&a*dPc_eq(i)/(de*Ctz*ga*ga*At)end do!计?算?各é燃?层?时骸?的?推?进?剂á燃?速ê值μrspeeddo i=0,num,1rspeed(i)=fun(Sa(i),Apa(i),key)*a*Pc_eq(i)**pnend do!计?算?各é燃?层?时骸?的?推?进?剂á的?推?力 值μFeqdo i=0,num,1Cf=0.0Cf=ga*sqrt(2*k*(1-(101325.0*Pe/Pc_eq(i))**((k-1)/k))/(k-1))Feq(i)=Cf*Pc_eq(i)*Atend do!计?算?燃?到?各é层?需â要癮的?时骸?间?值μtimetime=0.0do i=1,num,1time(i)=time(i-1)+de*2.0/(rspeed(i)+rspeed(i-1))end do!-----------------------------------------------------------!计?算?起e始?段?的?压1力 值μ、￠推?力 和æ时骸?间?!-----------------------------------------------------------do i=1,n,1P1_eq(i)=((Pc_eq(0)-Pig)/n)*(i-1)+Pigend dodo i=1,n,1Cf=0.0Cf=ga*sqrt(2*k*(1-(101325.0*Pe/P1_eq(i))**((k-1)/k))/(k-1))F1(i)=Cf*P1_eq(i)*Atend dodo i=1,n,1time1(i)=Apa(0)*len/(Ctz*ga*ga*At*(1-pn))*&log((Pc_eq(0)-Pig)/(Pc_eq(0)-P1_eq(i)))end do!-----------------------------------------------------------!计?算?加è上?起e始?段?时骸?工¤作痢?段?的?时骸?间?timedo i=0,num,1time(i)=time1(n)+time(i)end do!-----------------------------------------------------------!计?算?拖?尾2段?的?压1力 和æ时骸?间?值μ(按恪?照?绝?热â?膨ç胀æ过y程å计?算?,不?考?虑?有瓺剩骸?药?作痢?用?)!-----------------------------------------------------------do i=0,n-1,1P2_eq(i)=Pc_eq(num)-i*(Pc_eq(num)-101325.0)/(n-1) end dodo i=0,n-1,1Cf=0.0Cf=ga*sqrt(2*k*(1-(101325.0*Pe/P2_eq(i))**((k-1)/k))/(k-1)) F2(i)=Cf*P2_eq(i)*Atend dodo i=0,n-1,1time2(i)=time(num)+2.0*Apa(num)*len*&((P2_eq(i)/Pc_eq(num))**((1-k)/(2.0*k))-1.0)/&((k-1)*ga*At*sqrt(Roc*Toc))end do!将?平?衡a压1力 转羇化ˉ为a物?理え?大洙?气?压1的?形?式?do i=1,n,1P1_eq(i)=P1_eq(i)/101325.0end dodo i=0,num,1Pc_eq(i)=Pc_eq(i)/101325.0end dodo i=0,n-1,1P2_eq(i)=P2_eq(i)/101325.0end do!输?出?压1力 值μopen(10,file="export_inter_P.dat")do i=1,num+n+n-1,1if(i<=n-1) thenwrite(10,"(f15.5,2x,f15.5)") time1(i),P1_eq(i)else if(i>=n.and.i<num+n) thenwrite(10,"(f15.5,2x,f15.5)") time(i-n),Pc_eq(i-n) elsewrite(10,"(f15.5,2x,f15.5)") time2(i-num-n),P2_eq(i-num-n) end ifend doclose(10)!输?出?推?力 值μopen(20,file="export_inter_F.dat")do i=1,num+n+n-1,1if(i<=n-1) thenwrite(20,"(f15.5,2x,f15.5)") time1(i),F1(i)else if(i>=n.and.i<num+n) thenwrite(20,"(f15.5,2x,f15.5)") time(i-n),Feq(i-n)elsewrite(20,"(f15.5,2x,f15.5)") time2(i-num-n),F2(i-num-n) end ifend doclose(20)returnend subroutine!计?算?侵?蚀骸?比括?的?子哩?函ˉ数簓function fun(Sa,Apa,key)implicit noneinteger :: keyreal(kind=8) :: Sa,Apareal(kind=8) :: funif(key==1) thenif(Sa/Apa<=72.9) thenfun=1.0elsefun=1.3128-1.3249e-2*Sa/Apa+&1.5527e-4*(Sa/Apa)**2-4.3868e-7*(Sa/Apa)**3end ifelsefun=1.0end ifreturnend。