原因
所给链接的代码不完整，缺少ndd_fun函数。
顺便鄙视一下该链接的提供者，这么广为流传的东西下载居然还要积分，简直穷疯了。
代码
帮你好好找了一下，找到了完整的程序，供参考（全部代码保存到一个M文件运行即可，或直接下载附件）：
function ndd%59nian130
A=0.87 %q(炮)膛横断面积A dm^2
G=19%33.4 %d重 kg
W0=2.04 %药室容积 dm^3
l_g=25.0 %身管行程 dm
P_0 =30000 %起动压力 kpa
fai1=1.02 %次要功系数
K=1.03 %运动阻力系数φ1
theta =0.2 %火药热力系数
%=========================================
f=950000 %火药力 kg*dm/kg
alpha=1 %余容 dm^3/kg
delta=1.6 %火药重度γ
%==================================
ome=2.2%12.9 %第一种装药量森山 kg
u1=5.0024*10^-5 %第一种装药烧速系数 dm^3/(s*kg)
n1=0.82 %第一种装药的压力指数n1
lambda=-0.0071 %第一种装药形状特征量λ1
lambda_s=0 %第一种装药分裂点形状特征量λ1s
chi=1.00716 %第一种装药形状特征量χ1
chi_s=0 %第一种装药此举中分裂点形状特征量χ1s
mu=0 %第一种装药形状特征量μ1
et1=1.14*10^-2 %第一种装药药厚δ01
d1=2.5*10^-2 %第一种装药火药内径d1
Ro1=0 %药型系数α1
%=========================================
%常数与初值计算-----------------------------------------------------
l_0=W0/A
Delta=ome/W0
phi=K + ome/(3*G)
v_j=196*f*ome/(phi*theta*G)
v_j=sqrt(v_j)
B = 98*(et1*A)^2/( u1*u1*f*ome*phi*G )
B=B*(f*Delta)^(2-2*n1)
Z_s=1+Ro1*(d1/2+et1)/et1
p_0=P_0/(f*Delta)
psi_0=(1/Delta - 1/delta)/(f/P_0 + alpha - 1/delta)
Z_0=(sqrt(1+4*psi_0*lambda/chi) - 1)/(2*lambda)
%解算子------------------------------------------------------------
C = zeros(1,12)
C(1)=chiC(2)=lambdaC(3)=lambda_sC(4)=chi_sC(5)=Z_s%
C(6)=thetaC(7)=BC(8)=n1C(9)=DeltaC(10)=deltaC(11)=alphaC(12)=mu
C
y0=[Z_000psi_0]
options = odeset('outputfcn','odeplot')
[tt,y] = ode45(@ndd_fun,0:100,[Z_000],options,C)
l = y(:,2)
l = l*l_0
fl = find(l>=l_g)
fl = min(fl)
[tt,y] = ode45(@ndd_fun,0:0.005:fl,[Z_000],options,C)
Z = y(:,1)lx = y(:,2) vx = y(:,3)
psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...%%%%%%%%%
(Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) 答御) +...
(Z>=Z_s)*1
l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi
px = ( psi - vx.*vx )./( lx + l_psi )
p = px*f*Delta/100
v = vx*v_j/10
l = lx*l_0
t = tt*l_0*1000/v_j
fl = find(l>=l_g)
fl = min(fl)+1
p(fl:end)=[]v(fl:end)=[]l(fl:end)=[]t(fl:end)=[]
pd=px*f*Delta/100/(1+ome/3/fai1/G)
pt=pd*(1+ome/2/fai1/G)
aa=max(px)
M=find(px==aa)
Pm=[tt(M)*l_0*1000/v_j lx(M)*l_0 vx(M)*v_j/10 px(M)*f*Delta/100 pt(M) pd(M) psi(M) Z(M)]
%ll=length(tt)
ran=find(Z>=1)
ran=min(ran)
Zf=[tt(ran)*l_0*1000/v_j lx(ran)*l_0 vx(ran)*v_j/10 px(ran)*f*Delta/100 pt(ran) pd(ran) psi(ran) Z(ran)]
jie=find(psi>=1)
jie=min(jie)
psij=[tt(jie)*l_0*1000/v_j lx(jie)*l_0 vx(jie)*v_j/10 px(jie)*f*Delta/100 pt(jie) pd(jie) psi(jie) Z(jie)]
pg=[tt(end)*l_0*1000/v_j lx(end)*l_0 vx(end)*v_j/10 px(end)*f*Delta/100 pt(end) pd(end) psi(end) Z(end)]
Ry1=[ZfpsijpgPm]
Ry2=[tt*l_0*1000/v_j lx*l_0 vx*v_j/10 px*f*Delta/100 pt pd psi Z]
subplot(2,2,1)
plot(t,p,'linewidth',2)
grid on
xlabel('\fontsize{8}\bft (ms)')
ylabel('\fontsize{8}\bfp (kg/cm^{2})')
title('\fontsize{8}\bft-p曲线')
subplot(2,2,2)
plot(t,v,'linewidth',2)
grid on
xlabel('\fontsize{8}\bft (ms)')
ylabel('\fontsize{8}\bfv (m/s)')
title('\fontsize{8}\bft-v曲线')
subplot(2,2,3)
plot(l,p,'linewidth',2)
grid on
xlabel('\fontsize{8}\bfl (dm)')
ylabel('\fontsize{8}\bfp (kg/cm^{2})')
title('\fontsize{8}\bfl-p曲线')
subplot(2,2,4)
plot(l,v,'linewidth',2)
grid on
xlabel('\fontsize{8}\bfl (dm)')
ylabel('\fontsize{8}\bfv (m/s)')
title('\fontsize{8}\bfl-v曲线')
tspan = length(t)/20
tspan = 1:ceil(tspan):length(t)
tspan(end) = length(t)
fprintf(' t(ms) p(kg/cm^2) v(m/s) l(dm)')
format short g
Result = [t(tspan) p(tspan) v(tspan) l(tspan)]
format
%------------------------------------------------------------------
function dy = ndd_fun(t,y,C)
chi=C(1)lambda=C(2)lambda_s=C(3)chi_s=C(4)Z_s=C(5)mu=C(12)
theta=C(6)B=C(7)V=C(8)Delta=C(9)delta=C(10)alpha=C(11)
Z = y(1) l = y(2) v = y(3)
psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...
(Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...
(Z>=Z_s)*1
l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi
p = ( psi - v*v )/( l + l_psi )
dy(1) = sqrt(theta/(2*B))*(p^V)*(Z>=0&Z<=Z_s)
dy(2) = v
dy(3) = theta*p/2
dy = [dy(1)dy(2)dy(3)]
结果
d道d道有两种， 一是内部d道， 一是外部d道。
内d道讨论的是在d药击发后， d头离开q口前， 各种物理现象。 其中重要的有:
一. 膛压: 装药燃烧而扩张， 因为d头在前挡着， 机锁在后堵住， 在q膛中会产生极大的压力， 一般而言在数万磅/寸到数十万磅/寸之间。 这个压力是在d头脱离d茄罩知壳时， 推动d头的主要力量。 当然这个力量越大越好， 因为d头飞得越快， 在固定距离内， 受地心引力的影响越小。 但是如果不在适当的范围内， 也会产生发生危险。
二. 来复线: 来复线造成d头的旋转， 而使得d头的飞行稳定， 可是来复线的数量和线的快， 慢(快慢指的是来复线在多少长度完成360度旋转)和d头的重量有极重要的关系。 正确的d头用在适宜的来复线上， 会有较好的精确度。 例如说， 在使用 .223 的q里， 12 寸一圈以上的的， 适用55 gr 的d头。 如果是9寸一圈的， 就该用69 gr 以上的d头较好。
三. q管硬度: 在d头通过q管时， q管会像鞭子一样上下甩动。 动的幅度会影响到d头出口的位置。 同时， q管会发热， 金属因热而扩张， d头和来复线的密和度会受到影响。 要解决这个问题， 一般来说是增加q管的厚度。 因为增加厚度可以增加硬度而且减缓温度提高。
外d道主要是讨论d头出口后， 影响其飞行的各种因素。
任何在地球上的物体， 都会受到地心引力的影响。 (事实上光也会受到引力的影响， 但是光到底是波还是粒子， 还无定论)。 d头一出q口， 加速就停止了。 引力会将d头往地面拉。 所以任何d头的飞行路线都是弧形的。 如果q管与地面平行， d头永远不会和q管延长线的任何一点交会。 所以， q管都是微微朝上的。
d道与瞄准线示意图
这条弧线的弧度 (Trajectory)， 取决于d头出膛的初速和子d的流体系数 (co-efficient)。 初速大， d头在相等时间， 飞行距离远， 引力作用的时间短， 影响弧线的程度小， 飞行的弧线也就比较平坦。 平坦的d道表示d头不会偏离瞄准线太远， 对射击者而言， 简单的多了。 基本上是瞄那里就打那里， 不用担心调整准心或是调整瞄准点。
相信大家都听过归零这个名词， 主要就是因为d道和瞄准线的相对关系。 一个是直线， 一个是弧线， 最多只有两个交会点。 如何取决这两个交点， 就取决于q的用途了。 有光学瞄准镜的q， 如果有调整钮， 基本上是将外交点移动， 使d头准确的击中瞄准点。 军用步q， 因为在战斗中不会有时间去调整表尺。 交点多半是固定的。 例如 AK-47， 标准作业程序是将表尺设在两百公尺， 然后射击二十五公尺的标的。 这时可以调整表尺的上下左右， 当d头击在靶心时， 就完成了战斗归零。 此时这把q从零到三百公尺， 可以击中在瞄准点上下九寸之内。 对步兵而言， 这样的准确度已经够了。 想想看， 一个成年男子的上半身大约是十八寸， 宽约十寸。
舟尾形d头
d头闷带系数和d道的关系， 主要是因为d头必须在空气中飞行， 空气阻力会影响d头的飞行速度， 也就影响了d道。 一般来说， 流线形的d头有比较高的系数。 舟尾形 (boat-tail) 的d头比较不会受风力的影响而偏离瞄准线或是损失太多速度。 还有， 射击者往上和往下瞄准时， 都需要调整。 大约都是要瞄高一些。 这里有个简单的公式。
30 度角 (朝上或是朝下都是一样):
距离X 0.7
45 度角
距离X 0.9
举例来说， 600 码的山坡上有一个目标。 射手必须朝上仰角三十度射击， 经过计算， 600X0.7=420， 射手要调整表尺或是瞄准镜， 有如要射击一个420 码的目标一样。 这只是一个估计值， 不是绝对精确的。
有兴趣的话， 可以买一台测速机 (Chronograph)， 测量d头的出口速度， 另外还有一些d道软体， 可以绘出d道的理论值颤消。 是很好的参考资料， 可以帮助射击的精确度
1、起始d道，研究从击发开始到d丸的d带全部挤进膛线，或克服其他起动阻力这一阶段中，点火药的返兆庆点火与传火、装药燃烧、d丸挤进过程以及压力漏握波的形成与发展等规律及有关现象。
2、膛内d道，研究从d丸的d带全部挤进膛线到d丸飞出膛口这一阶段的d丸运动、火药燃烧、物质流动以及能量转换等规律及有关现象。
3、中间d道，研究从d丸出膛口到脱离火药燃气的力学影响这一阶段的膛口流场对d丸运动规律的影响，以及伴随膛内火药燃气排空过程发生的有关现象。
4、膛外d道，研究d丸在脱离膛口流场影响猜颂之后，在空中飞行的运动规律以及有关现象。
5、终点d道，研究d丸在目标区域发生的现象与运动规律，对目标的作用(如爆炸、冲击、侵彻等)机理及威力效应等。
扩展资料：
d道导d在火箭发动机推力作用下按预定程序飞行，关机后按自由抛物体轨迹飞行的导d。其飞行d道一般分为主动段和被动段。
主动段（又称动力飞行段或助推段）是导d在火箭发动机推力和制导系统作用下，从发射点起飞到火箭发动机关机时的飞行路径；
被动段包括自由飞行段和再入段，是导d按照在主动段终点获得的给定速度和d道仪角作惯性飞行，到d头起爆的路径。
参考资料：百度百科-d道
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）