有限差分法计算固体火箭发动机内弹道

!计?算?星?型í装痢?药?的?几?何?尺?寸?!参?数簓符?号?说μ明÷!-----------------------------------------------------------------------------------!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)-foresmaxif(error>=0.0001)thenthet=thet+0.00001else if(error<=-0.0001)thenthet=thet-0.00001elseexitend 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**2write(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)!调獭?用?子哩?程ì序ò计?算?装痢?药?内ú弹獭?道台?曲ú线?callinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp) !输?出?星?型í装痢?药?的?几?何?参?数簓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!该?子哩?程ì序ò用?于 ?计?算?零?维?变?截?面?燃?烧?装痢?药?的?内ú弹獭?道台?!利?用?有瓺限T差?分?法ぁ?计?算?内ú弹獭?道台?曲ú线?!可é用?于 ?计?算?侵?蚀骸?燃?烧?效§应畖下?的?内ú弹獭?道台?曲ú线?!--------------------------------------------------------------------------!d外猘径?,len长¤度è,e1平?均ù肉╝厚?,n星?角?数簓,thet星?边?夹D角?,epsilon角?度è系μ数簓!r过y度è圆2弧?半?径?,r1星?角?圆2弧?半?径?和íl特?征÷尺?寸?,key表括?示?是?否?考?虑?侵?蚀骸?燃?烧?!density_p推?进?剂á的?密ü度è,k比括?热è?比括?ga系μ数簓,C特?征÷速ù度è,mpeff药?柱ù质ê量?!a速ù度è系μ数簓,pn压1强?指?数簓,At喉í部?面?积y,Pc燃?烧?室酣?的?设Θ?计?压1力 ,Isp 表括?示?理え?论?比括?冲?!dt时骸?间?步?长¤,e燃?层?厚?度è,time时骸?间?,Sa同?一?时骸?刻ì轴á向ò的?各÷节ú点?的?燃?面?,Apa通 ?气?面?积y!Poc各÷节ú点?的?压1力 ,Pe喷?管ü出?口ú压1力 ,P_I压1力 冲?量?,ep侵?蚀骸?比括?!P_av平?均ù压1力 ,F_av平?均ù推?力 ,I0总哩?冲?,Im重?量?比括?冲?Iv体?积y比括?冲?,Cf 推?力 系μ数簓!--------------------------------------------------------------------------subroutineinternal_ballistics0(d,len,e1,n,thet,epsilon,r,r1,l,a,pn,Ctz,At,k,density_p,mpeff,Isp) implicit nonereal,parameter::Pi=3.14integer::n,keyreal(kind=8),intent(in)::d,len,e1,thet,epsilon,r,r1,lreal(kind=8),intent(in)::density_p,k,Ctz,mpeff,Ispreal(kind=8)::a,pn,At,dt,e,time,Sa,Apa,Poc,Pe,P_I,Freal(kind=8)::ep,P_av,F_av,I0,Im,Iv,Cf,gaPe=101325.0P_I=0.0!药?柱ù的?能ü量?特?性?参?量?读á入?open(3,file="inter_parameter.dat")read(3,*)keyread(3,*)dtread(3,*)Pocclose(3)!计?算?点?火e压1强?时骸?的?推?力 值μPoc=Poc*101325.0ga=sqrt(k)*(2/(k+1))**((k+1)/(k-1)/2)Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Ate=0.0time=0.0open(10,file="export_Poc.dat")open(20,file="export_F.dat")write(10,"(f8.5,2X,f10.5)")time,Poc/101325.0write(20,"(f8.5,2X,f15.5)")time,Fdo while(.true.)!计?算?燃?面?和í通 ?气?面?积yif(e<=r1)thenSa=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-e)*Pi/n/l)*l*len Apa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-&cotan(thet/2))/l**2+n*(r1-e)**2*(thet/2+cotan(thet/2)-Pi/2)/l**2)*l**2else if(e>r1.and.e<=(l*sin(epsilon*Pi/n)/cos(thet/2)-r))thenSa=2*n*(sin(epsilon*Pi/n)/sin(thet/2)+(1-epsilon)*Pi/n+&(e+r)*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l)*l*lenApa=(n*((1-epsilon)*Pi/n+sin(epsilon*Pi/n)*(cos(epsilon*Pi/n)-&sin(epsilon*Pi/n)*cotan(thet/2)))+2*n*(r+e)*(sin(epsilon*Pi/n)/&sin(thet/2)+(1-epsilon)*Pi/n)/l+n*(r+e)**2*(Pi/n+Pi/2-thet/2-cotan(thet/2))/l**2)*l**2 else if(e>(l*sin(epsilon*Pi/n)/cos(thet/2)-r).and.e<=e1)thenSa=2*n*((1-epsilon)*Pi/n+(r+e)*(Pi/n+asin(l*sin(epsilon*Pi/n)/(e+r)))/l)*l*len Apa=n*((1-epsilon)*Pi*(1+(r+e)/l)**2/n+sin(epsilon*Pi/n)*(sqrt((r+e)**2/l**2-& sin(epsilon*Pi/n)**2)+cos(epsilon*Pi/n))+(r+e)**2*&(epsilon*Pi/n+asin(l*sin(epsilon*Pi/n)/(r+e)))/l**2)*l**2elseSa=0.0end if!计?算?侵?蚀骸?比括?if(key==1)thenif(Sa/Apa<=72.9)thenep=1.0elseep=1.3128-1.3249e-2*Sa/Apa+&1.5527e-4*(Sa/Apa)**2-4.3868e-7*(Sa/Apa)**3end ifelseep=1.0end if!采é用?有瓺限T差?分?法ぁ?计?算?压1强?Poc=Poc+(ga*ga*Ctz*Ctz*dt/(Apa*len))*&(density_p*Sa*ep*a*Poc**pn-Poc*At/Ctz)Cf=ga*sqrt(2*k*(1-(Pe/Poc)**((k-1)/k))/(k-1))F=Cf*Poc*Attime=time+dte=e+dt*ep*a*Poc**pnif(e>e1.and.Poc<101325.0)exitwrite(10,"(f8.5,2X,f10.5)")time,Poc/101325.0write(20,"(f8.5,2X,f15.5)")time,FP_I=P_I+Poc*dtend doclose(10)close(20)!循-环·结á束?P_av=P_I/timeCf=ga*sqrt(2*k*(1-(Pe/P_av)**((k-1)/k))/(k-1))F_av=Cf*P_av*AtI0=F_av*timeIm=I0/mpeff/9.81Iv=I0/(Apa*len)!输?出?固ì体?火e箭y发ぁ?动ˉ机ú的?工¤作痢?参?数簓值μopen(30,file="export_propellant.dat")write(30,"('平?均ù压1力 p_av=',f9.5)")p_av/101325.0write(30,"('平?均ù推?力 F_av=',f10.5,1X,'KN')")F_av/1000.0 write(30,"('总哩?冲?I0=',f10.5,1X,'KN.s')")I0/1000.0write(30,"('重?量?比括?冲?Im=',f10.5,1X,'s')")Imwrite(30,"('体?积y比括?冲?Iv=',f10.5,1X,'KN.s/m3')")Iv/1000.0 write(30,"('工¤作痢?时骸?间?time=',f8.5,1X,'s')")timeclose(30)returnend subroutine。