内弹道计算
- 格式:doc
- 大小:69.00 KB
- 文档页数:5
内弹道是指射程较短的导弹或火箭弹在飞行过程中受到大气阻力和重力等作用的飞行轨迹。
内弹道理论研究的是导弹或火箭弹在发射后到离开大气层再进入大气层末时的飞行过程。
内弹道包括导弹或火箭弹在发射后的加速、稳定、制导、飞行以及飞行过程中的动力学性能仿真等诸多内容。
内弹道有着复杂的飞行特性和动力学方程,在实际工程中需要进行准确的计算和仿真。
内弹道的计算中,龙格库塔(Runge-Kutta)法是一种常用的数值积分方法,在求解微分方程等领域有着广泛的应用。
龙格库塔法是由数学家奥特翁格(C. W. Runge)和马丁庫塔(M. W. J. Kutta)于1900年提出的,用于求解常微分方程初值问题,其优点是精度较高,适用范围广。
在内弹道计算中,可以利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算,得到较为准确的飞行轨迹数据。
在实际工程中,为了方便进行内弹道的计算,可以使用Matlab等数学建模和仿真软件。
Matlab是一种常用的科学计算软件,具有强大的数值计算和仿真功能,可以用于内弹道计算中的龙格库塔法数值模拟。
在Matlab中,可以编写相应的程序,利用龙格库塔法对导弹或火箭弹的飞行过程进行仿真和计算,得到准确的飞行轨迹和动力学性能数据。
内弹道计算是导弹或火箭弹研究设计中的重要内容,龙格库塔法是一种常用的数值积分方法,Matlab是一种常用的科学计算软件,它们的应用能够有效地进行内弹道的计算和仿真,为导弹或火箭弹的研制提供重要的技术支持。
随着技术的不断发展,内弹道计算已经成为导弹或火箭弹研究设计中不可或缺的一部分。
在内弹道计算中,龙格库塔法是一种常用的数值积分方法,可以对导弹或火箭弹的飞行轨迹进行数值模拟和计算,提供准确的飞行轨迹数据。
而Matlab作为一种强大的科学计算软件,对于内弹道的计算和仿真也有着重要的应用价值。
在实际工程中,使用Matlab编写程序,利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算,将为导弹或火箭弹的研制提供重要的技术支持。
火炮内弹道求解与计算
火炮内弹道是指火炮射击时炮弹在火炮内的运动轨迹。
要解决火炮内弹道问题,需要考虑炮弹在炮管内的运动特性,以及发射药燃烧产生的气体对炮弹的推动力。
本文将从炮弹的运动方程入手,分析火炮内弹道的解法并进行计算。
炮弹的运动方程可以表示为:
ma = F - mg - fd - fL
其中m是炮弹的质量,a是炮弹在炮管内的加速度,F是发射药燃烧产生的推动力,g是重力加速度,fd是炮弹在炮管内受到的阻力,fL是炮弹在炮管内受到的气体偏转力。
在火炮运动方程中,炮弹在炮管内的加速度a是常量,可以通过测量炮弹的初速度和射程得到。
炮弹的初速度可以通过实验或者计算得到。
发射药燃烧产生的推动力F可以通过推进药的燃烧速率和燃烧产物的排放速度进行计算。
通过实验或者模拟可以得到推进药的燃烧速率和燃烧产物的排放速度。
炮弹在炮管内受到的阻力fd可以通过火炮内管壁的摩擦力和火药燃烧产生的气体对炮弹的阻力进行计算。
火炮内管壁的摩擦力可以由实验和数学模型得到。
火药燃烧产生的气体对炮弹的阻力可以通过实验和气体动力学模型计算。
炮弹在炮管内受到的气体偏转力fL可以通过气体对炮弹的作用力和炮弹的偏转角度进行计算。
气体对炮弹的作用力可以由实验和气体动力学模型得到。
炮弹的偏转角度可以由实验或者数学模型计算。
通过解决火炮内弹道问题,可以得到炮弹的运动轨迹和射程。
在实际应用中,可以通过对火炮内弹道进行数值模拟和优化计算,提高火炮的射击精度和射程。
舰载蒸汽弹射内弹道设计计算舰载蒸汽弹射内弹道设计计算舰载蒸汽弹射是现代航母起飞的最常用方式之一,它通过利用高压蒸汽推动喷气式飞机飞出航母甲板,具有快速高效和适应各种飞机的特点。
内弹道设计计算是舰载蒸汽弹射系统设计的重要部分,通过准确计算飞机的起飞质量、速度和加速度等参数,以及考虑飞行姿态和气动特性,从而确保安全、稳定和高效的起飞过程。
一、舰载蒸汽弹射系统工作原理舰载蒸汽弹射系统是由蒸汽动力机组、蒸汽管路、弹射准备、准备、发射准备控制系统等组成的。
飞机进入弹射器后,与弹射器碰触的瞬间,弹射器向后推出一进气孔以外的压缩空气,压缩空气进入涡轮机发生回转作用。
二、内弹道设计计算1. 起飞重量计算起飞重量是指飞机在起飞时的总重量,包括机身、燃料、弹药、载荷和人员等。
起飞重量的计算是内弹道设计计算的重要基础。
其计算公式如下:起飞重量 = 机身重量 + 最大燃油重量 + 载荷 + 弹药 + 人员2. 加速度计算加速度是弹射过程中比较关键的参数,其大小直接决定飞机的起飞速度和高度。
其计算公式如下:加速度 = 2 * 起飞总推力 / 起飞重量起飞总推力包括飞机引擎产生的推力和蒸汽弹射系统提供的推力。
一般情况下,弹射器的起飞总推力要达到飞机重量的1.2倍以上,以确保飞机在起飞过程中有足够的加速度。
3. 起飞速度计算起飞速度是指飞机在弹射器上达到准备起飞状态所需的速度,取决于加速度、飞机重量和气动特性等因素。
根据实际情况,起飞速度一般在200至250节之间。
其计算公式如下:起飞速度= √(2 * 起飞重量 * 加速度 / 飞机空气阻力系数 * 高度密度)飞机空气阻力系数和高度密度是通过实验和理论计算得出的参数。
4. 起飞高度计算起飞高度是指飞机在离开航母甲板时的高度,并直接关系到飞机在起飞过程中的安全和稳定。
其计算公式如下:起飞高度 = 起飞速度 * 弹射器长度弹射器长度是通过实际测量得出的参数,通常在80至100米之间。
59-130加农炮内弹道计算function ndd%59-130A=1.394; %枪(炮)膛横断面积A dm^2G=33.4; %弹重kgW0=18.56; %药室容积dm^3l_g=59.52; %身管行程dmP_0 =30000; %起动压力kpafai1=1.02; %次要功系数K=1.03; %运动阻力系数φ1theta =0.2; %火药热力系数%=========================================f=950000; %火药力kg*dm/kgalpha=1; %余容dm^3/kgdelta=1.6; %火药重度γ%==================================ome=12.9; %第一种装药量kgu1=5.0024*10^-5; %第一种装药烧速系数dm^3/(s*kg)n1=0.82; %第一种装药的压力指数n1lambda=-0.0071; %第一种装药形状特征量λ1lambda_s=0; %第一种装药分裂点形状特征量λ1schi=1.00716; %第一种装药形状特征量χ1chi_s=0; %第一种装药分裂点形状特征量χ1smu=0; %第一种装药形状特征量μ1et1=1.14*10^-2; %第一种装药药厚δ01d1=2.5*10^-2; %第一种装药火药内径d1Ro1=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;% ********************* ndd- fun*********************** 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)];二.运行结果Result =t(ms) p(kg/cm^2) v(m/s) l(dm)0 300 0 00.65 500.1 9.022 0.026841.3 792.68 23.657 0.129411.95 1192.1 46.219 0.351562.6 1690.7 79.13 0.752713.25 2242 124.22 1.40683.9 2759.5 181.76 2.39474.55 3146.1 249.84 3.79255.2 3343.5 324.65 5.65725.85 3356.6 401.84 8.01816.5 3233.7 477.64 10.8787.15 3033.1 549.63 14.2197.8 2801.1 616.59 18.0128.4499 2566.8 678.16 22.2239.0999 2345.5 734.48 26.8179.7499 1968.6 783.83 31.75610.4 1671.4 825.45 36.9911.05 1437.9 861.01 42.47411.7 1251.7 891.79 48.17212.939 988.32 940.43 59.54。
大口径机枪内弹道及导气室压力的计算1. 前言大口径机枪是当今最先进的军事武器之一,其弹道精度和射程都在不断提高。
在机枪的设计和制造过程中,数学模型的建立和弹道计算是非常重要的一环。
本文将探讨大口径机枪内弹道及导气室压力的计算方法。
2. 大口径机枪内弹道的计算大口径机枪内弹道的计算需要考虑多种因素,包括弹丸的质量、速度、空气阻力、自旋等。
这些因素都会对弹道轨迹产生影响,因此需要建立数学模型,通过模拟计算得到准确的弹道轨迹。
2.1 弹道轨迹计算模型弹道轨迹计算模型分为两类:解析方法和数值方法。
解析方法是通过数学公式推导得出弹道轨迹方程,可以得到准确的解析解,但适用范围有限,只适用于某些特定情况。
数值方法是通过离散化求解弹道轨迹,在计算机上进行数值仿真,适用范围更广,但是需要计算机进行复杂的运算,速度较慢。
对于大口径机枪内弹道的计算,常用的数值方法是欧拉法和龙格-库塔法。
欧拉法是通过离散化时间对连续的弹道轨迹进行拟合,在每个时间步长内计算弹丸的速度和位置。
龙格-库塔法是欧拉法的高阶拓展,通过更复杂的公式计算弹丸的速度和位置,精度更高。
2.2 弹道轨迹计算流程弹道轨迹计算的流程包括以下步骤:1. 建立弹道轨迹模型,确定算法、时间步长、精度。
2. 输入弹丸初速度、大小、自旋、质量等参数。
3. 计算弹丸的速度和位置,考虑空气阻力、引力等因素。
4. 判断弹丸是否击中目标,如果没有继续计算,直到弹丸击中目标或离开射程。
5. 输出计算结果,包括弹道轨迹、飞行时间、命中点等信息。
2.3 建模与计算工具建模与计算工具包括Mathematica、Matlab、Python等,其中Python是最常用的工具之一,因为Python语言简洁易学,有众多科学计算库和开源工具包。
3. 导气室压力的计算导气室是大口径机枪内的一个重要组成部分,可以提高弹丸的初速度和射程。
导气室压力的计算是确定导气室的设计和参数的关键步骤。
3.1 导气室原理导气室是由枪管和枪膛组成的,当枪管中的火药燃烧产生高温高压气体时,这些气体进入导气室,使其内压力急剧增加,从而产生推动力,使弹丸加速并射出。
后抛式弹丸空中抛射内弹道计算引言:抛式弹丸空中抛射内弹道计算是研究抛式武器弹丸在空中抛射过程中的弹道轨迹的一项重要内容。
通过对抛射角度、初速度、重力等因素的计算分析,可以确定弹丸的轨迹,为军事和工程领域的设计和应用提供理论依据。
一、抛式弹丸的基本概念抛式弹丸是指通过一定的装置将弹丸以一定的角度和初速度抛射到空中的武器系统。
在这个过程中,弹丸受到重力的作用,同时还可能受到空气阻力等因素的影响。
内弹道是指弹丸在抛射后离开抛射装置进入空中,到达弹道稳定状态之前的过程。
二、抛式弹丸空中抛射内弹道计算的基本原理抛式弹丸空中抛射内弹道计算主要基于牛顿力学和空气动力学的基本原理。
在计算过程中,需要考虑以下几个因素:1. 抛射角度:抛射角度是指弹丸的抛射方向与水平方向之间的夹角。
不同的抛射角度会影响弹丸的最大射程和轨迹形状。
2. 初速度:初速度是指弹丸离开抛射装置时的速度。
初速度的大小与弹丸的射程和轨迹形状密切相关。
3. 重力:重力是指地球对弹丸的作用力。
重力的大小与地球的重力加速度和弹丸的质量有关。
4. 空气阻力:空气阻力是指弹丸在运动过程中,由于与空气分子的碰撞产生的阻力。
空气阻力的大小与弹丸的速度、形状和密度等因素有关。
5. 其他因素:还有一些其他因素也会影响弹丸的轨迹,如风速、风向、弹丸的自旋等。
三、抛式弹丸空中抛射内弹道计算的方法抛式弹丸空中抛射内弹道计算可以使用数值计算方法或解析计算方法。
数值计算方法通常利用计算机模拟弹丸的运动过程,通过迭代计算得到弹丸的位置和速度。
解析计算方法则通过建立数学模型,利用数学方程求解弹丸的轨迹方程。
1. 数值计算方法:数值计算方法主要包括欧拉法、龙格-库塔法等。
这些方法通过将抛式弹丸的运动过程离散化,将连续的运动变为离散的运动,然后利用迭代算法计算弹丸在每个时间步长上的位置和速度,从而得到整个弹丸的轨迹。
2. 解析计算方法:解析计算方法主要通过建立数学模型,利用微分方程或积分方程求解弹丸的轨迹方程。
发射动力系统内弹道优化设计计算发射动力系统内弹道优化设计计算发射动力系统内弹道优化设计计算是探索重点任务之一,因为它关系到弹道导航与控制系统的精度和可靠性,直接影响到导弹的打击效果和命中率。
本文将对发射动力系统内弹道优化设计计算进行详细介绍。
一、发射动力系统内弹道发射动力系统内弹道是导弹在离开发射台后到达目标点之前的轨迹。
一般来说,内弹道采用了三段加速法,即在离发射台距离较远的位置采用第一段加速,使导弹进入空气稀薄层中加速追踪目标;在距离目标点较远的位置采用第二段加速;在离目标点较近的位置采用第三段加速,使导弹能够击中目标,实现任意角度的攻击。
二、内弹道优化设计计算内弹道的优化设计计算目的是确定最佳的飞行计划和调整飞行参数,以使导弹能够以最小的时间、最小的燃料消耗和最大的精度击中目标。
(一)导引律选择导引律是导弹内弹道控制系统的核心,选择合适的导引律可以有效提高导弹的命中率和抗干扰性能。
常见的导引律有比例导引律、比例修正导引律、比例-积分导引律和预测导引律等。
在具体设计时需要根据目标类型、干扰环境和系统要求等综合因素进行选择。
(二)控制极点设计内弹道控制极点的设计是使导弹飞行稳定、准确的保证,控制极点对内弹道的稳定性、敏感度和过冲量等指标起到直接的影响。
调节控制极点的位置和数量可以精确控制导弹的动态行为,如响应速度、阻尼比、稳定性和过冲量等参数。
(三)预测法控制预测法控制是一种高级的弹道控制方法,与常规的比例-积分导引律不同的是,它使用预测技术来基于中间目标预测趋势,根据预测结果对导弹控制系统进行修正,使导弹能够更快、更准确地找到目标。
预测法控制可以提高导弹的抗干扰能力和命中率,特别适用于高速飞行和大气干扰条件下的导弹控制。
(四)弹体设计弹体设计是导弹内弹道优化设计的重要环节,它涉及到空气动力学、力学和材料科学等多学科交叉领域。
弹体设计的关键在于降低弹体的阻力和重量,提高弹体的机动性和抗干扰性能。
59-130加农炮内弹道计算function ndd%59-130A=1.394; %枪(炮)膛横断面积A dm^2G=33.4; %弹重kgW0=18.56; %药室容积dm^3l_g=59.52; %身管行程dmP_0 =30000; %起动压力kpafai1=1.02; %次要功系数K=1.03; %运动阻力系数φ1theta =0.2; %火药热力系数%=========================================f=950000; %火药力kg*dm/kgalpha=1; %余容dm^3/kgdelta=1.6; %火药重度γ%==================================ome=12.9; %第一种装药量kgu1=5.0024*10^-5; %第一种装药烧速系数dm^3/(s*kg)n1=0.82; %第一种装药的压力指数n1lambda=-0.0071; %第一种装药形状特征量λ1lambda_s=0; %第一种装药分裂点形状特征量λ1schi=1.00716; %第一种装药形状特征量χ1chi_s=0; %第一种装药分裂点形状特征量χ1smu=0; %第一种装药形状特征量μ1et1=1.14*10^-2; %第一种装药药厚δ01d1=2.5*10^-2; %第一种装药火药内径d1Ro1=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;% ********************* ndd- fun*********************** 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)];二.运行结果Result =t(ms) p(kg/cm^2) v(m/s) l(dm)0 300 0 00.65 500.1 9.022 0.026841.3 792.68 23.657 0.129411.95 1192.1 46.219 0.351562.6 1690.7 79.13 0.752713.25 2242 124.22 1.40683.9 2759.5 181.76 2.39474.55 3146.1 249.84 3.79255.2 3343.5 324.65 5.65725.85 3356.6 401.84 8.01816.5 3233.7 477.64 10.8787.15 3033.1 549.63 14.2197.8 2801.1 616.59 18.0128.4499 2566.8 678.16 22.2239.0999 2345.5 734.48 26.8179.7499 1968.6 783.83 31.75610.4 1671.4 825.45 36.9911.05 1437.9 861.01 42.47411.7 1251.7 891.79 48.17212.939 988.32 940.43 59.54。