Question

matlab 内弹道计算问题

如题，有matlab计算内弹道计算主程序，但运行时出错，新人求解，报错在@ndd_fun等三处地方，主程序地址http://wenku.baidu.com/link?url=qEcPqaO6l0Q9W1dnhTJys0nrK9qbwQ9LhvE0gTRi80f1LDlUxnQdW0tVwQyURvpWYQaf-VkVyFF5ZF_bIyEjh9xPzMeNXsJtbl0_RaGzqqO

Answer 1

原因

所给链接的代码不完整，缺少ndd_fun函数。

顺便鄙视一下该链接的提供者，这么广为流传的东西下载居然还要积分，简直穷疯了。

 

代码

帮你好好找了一下，找到了完整的程序，供参考（全部代码保存到一个M文件运行即可，或直接下载附件）：

function ndd
%59nian130
A=0.87;         %枪(炮)膛横断面积A  dm^2
G=19;%33.4;           %弹重  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)=chi;C(2)=lambda;C(3)=lambda_s;C(4)=chi_s;C(5)=Z_s;%
C(6)=theta;C(7)=B;C(8)=n1;C(9)=Delta;C(10)=delta;C(11)=alpha;C(12)=mu;
C;
y0=[Z_0;0;0;psi_0];
options = odeset('outputfcn','odeplot');
[tt,y] = ode45(@ndd_fun,0:100,[Z_0;0;0],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_0;0;0],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=[Zf;psij;pg;Pm];
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)];

 

结果

 

Answer 2

那个代码里面没有给ndd_fun这个函数.
这个主要是常微分方程组的表达式.
涉及到专业的东西,真不好说这个函数的具体内容.

你最好联系一下原作者.他代码没有给全.