固体火箭发动机零维内弹道的数值计算

        在上火发课程的时候，老师提了一嘴零维内弹道的基本方程没法用龙格库塔法，需要做近似或者用软件模拟。头铁的我想了很长时间，终于想到了用龙格库塔法硬解这个微分方程的手段。在这篇博客记录一下我是如何完成固体火箭发动机零维内弹道计算的大作业，并附上了完整的代码，也许能帮助学弟学妹少走些弯路，快些弄完这门课的大作业。

        固体火箭发动机零维内弹道的微分方程为：

         其中，V和Ab都是关于e的函数。而e和t的关系为：

         很显然，如果将V和Ab用t表示，上述微分方程中就出现了积分。这当然无法使用龙格库塔了。因此需要做一些变形。

        不妨对p对t的导数下手：

        这样，微分方程中就只剩下p,e两个未知数了，而且没有积分项，因此可以使用龙格库塔法。

        微分方程变形为：

         好了，这个形式完全能使用龙格库塔法。

        在利用龙格库塔法求出积分后，我们知道了p=p(e)。但我们最终的目标时求出p=p(t)。此时的p=p(t)依旧无法求出，但我们已经可以求取t=t(e)了。

        由于：

         这里就可以使用复化梯形来计算t=t(e)了。

        在求出t=t(e)和p=p(e)之后，我们只需取不同的e值，分别计算对应的p(e)和t(e),就可以得到p(t)了。

        得到的p=p(t)如下：

        完整代码如下：

#coding:gb2312
import numpy as np
import matplotlib.pyplot as plt

######这里开始是定义一些常数，便于写后续代码
T_k=0.6528
CRoA=1.959*10**7
CAt=2.7325
ed=27.25*10**(-3)*(np.sin(0.8*np.pi/6))/(np.cos(32*np.pi/180))-2.8*10**(-3)
k1=(0.02725*np.sin(0.8*np.pi/6))/(np.sin(32*np.pi/180))
k2=0.02725*0.2*np.pi/6
k3=np.pi/2+np.pi/6-32*np.pi/180-1/np.tan(32*np.pi/180)
k4=0.02725*np.pi/6*0.2
k5=0.5*(0.02725**2)*0.2*np.pi/6 
k6=0.5*(np.cos(0.8*np.pi/6)-np.sin(0.8*np.pi/6)/np.tan(32*np.pi/180))*(0.02725**2)*np.sin(0.8*np.pi/6)
k7=np.sin(0.8*np.pi/6)/np.sin(32*np.pi/180)+0.2*np.pi/6
k8=np.pi/2+np.pi/6-32*np.pi/180-1/np.tan(32*np.pi/180)
k9=0.02725**2/2
k10=0.2*np.pi/6
k11=np.sin(0.8*np.pi/6)
k12=np.cos(0.8*np.pi/6)
k13=0.8*np.pi/6
r=0.0028
l=0.02725
a=1.12356*10**-5
R=0.109/2
ypsl=0.8

#该函数的功能是求燃面
def Ab(e,type):
    num=e.shape[0]
    arr=np.zeros(num)
    if(type==0):
        arr=(k1+k2+(e+r)*k3)
    elif(type==1):
        arr=(k4+(e+r)*(np.pi/6+np.arcsin((0.02725/(e+r))*np.sin(0.8*np.pi/6))))
    else:
        arr=(e+r)*(np.arccos(((e+r)**2+l**2-R**2)/(2*(e+r)*l))-np.pi/2+ypsl*np.pi/6-np.arccos(l/(e+r)*np.sin(ypsl*np.pi/6)))
    return arr*12*1.1

#该函数的功能是求V
def V(e,type):
    num=e.shape[0]
    arr=np.zeros(num)
    if(type==0):
        arr=12*(k5+k6+0.02725*(e+r)*k7+0.02725/2*((e+r)**2)*k8)
    elif(type==1):
        arr=(12*k9*((1+(e+r)/l)**2*k10+k11*(np.sqrt(((e+r)/l)**2-k11**2)+k12)+((e+r)/l)**2*(k13+np.arcsin(k11/((e+r)/l)))))
    else:
        arr=e*0+0.00891
    return arr

########这里开始的十多行是预先分配好燃烧掉肉厚/燃面/体积 三个变量对应数组的内存
elen=100000
e=np.zeros(3*elen)
ab=np.zeros(3*elen)
v=np.zeros(3*elen)
e[0:elen]=np.linspace(0,0.0102695,elen)
ab[0:elen]=Ab(e[0:elen],0)
v[0:elen]=V(e[0:elen],0)
e[elen:2*elen]=np.linspace(0.01026953,0.02445,elen)
ab[elen:2*elen]=Ab(e[elen:2*elen],1)
v[elen:2*elen]=V(e[elen:2*elen],1)
e[elen*2:3*elen]=np.linspace(0.02445,0.02881,elen)
ab[elen*2:3*elen]=Ab(e[elen*2:3*elen],2)
v[elen*2:3*elen]=V(e[elen*2:3*elen],2)

#这个函数是等式 dy/de=……一式中右侧的部分，便于后续使用龙格库塔法
def f(nume,p):
    a1=(1.78632*10**9*ab[nume])/(v[nume])
    b1=(103629*p**0.5635)/(v[nume])
    return a1-b1

#这几个h是龙格库塔法的步长。在不同阶段,步长是不相同的。
h1=2.0*ed/elen
h2=2.0*(0.02445-ed)/elen
h3=2.0*(0.02881-0.02445)/elen
y=[300000]
E=[e[0]]

#这个for循环是利用四阶龙格库塔法求出p=p(e)
for i in range(int(elen*3/2)-1):
    if(i>elen/2):
        h1=h2
    if(i>elen):
        h1=h3
    E.append(e[2*i+2])
    t1=f(2*i,y[i])
    t2=f(2*i+1,y[i]+h1/2*t1)
    t3=f(2*i+1,y[i]+h1/2*t2)
    t4=f(2*i+2,y[i]+h1*t3)
    yn_1=y[i]+h1/6*(t1+2*t2+2*t3+t4)
    y.append(yn_1)
#接下来的几行是准备好r数组 
y=np.array(y)
y=y
r=a*y**0.4365
dt_de=1/r
t=[0]
i=0
#这个循环是利用复化梯形算法来求出t=t(e)
while i<len(E)-1:
    
    temp=t[i]+(E[i+1]-E[i])*(dt_de[i]+dt_de[i+1])/2
    t.append(temp)
    i=i+1
#剩下的是画图
plt.plot(E,y/1000000)
plt.xlabel('e/m')
plt.ylabel('pressure/MPa')

plt.show()

plt.plot(E,t)
plt.xlabel('e/m')
plt.ylabel('time/s')
plt.show()

plt.plot(t,y/1000000)
plt.xlim(-0.1,3.4)
plt.xlabel('time/s')
plt.ylabel('pressure/MPa')
plt.show()