| 查看: 1062 | 回复: 6 | |||
[求助]
求大神帮忙看一下程序，为什么一直无法执行选择之后的循环 已有2人参与
|
|
编了一个程序如下，希望p<pjc时使用一个循环，p>=pjc时使用另一个循环。在编程过程中发现不管是使用if选择语句还是while循环语句都是只能执行p<pjc时的循环，当p>pjc时程序就停止了，请各位大神帮忙看一下是什么原因，谢谢大家。 clear all p0=101.3;kp=7.9;T0=300;md0=1.3;u0=0.5;kz=4; Tbv=2221;Tbp=1867; V=5000;kl=1.4;kr=1.25; pb=p0*kp; pjc=p0+5; t(1)=0;tj=0.001; Vr(1)=0; Vrj(1)=0; p(1)=p0; u(1)=u0; krj(1)=0; rr(1)=0; mdw(1)=md0; Tw(1)=T0; i=1; %开始循环计算 while p(i)<=pb & p(i)>=p0 while p<=pjc & i<319 i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1); krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl)); p(i)=p(i-1)+krj(i)*(pb-p0); Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl)); Vw(i)=V-Vr(i); rr(i)=((3*Vr(i))/(4*pi))^(1/3); u(i)=((rr(i)-rr(i-1))/tj); mdw(i)=md0*((p(i)/p0)^(1/kl)); Tw(i)=T0*((p(i)/p0)^((kl-1)/kl)); Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i)); end while p>pjc i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1); krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl)); p(i)=p(i-1)+krj(i)*(pb-p0); pb(i)=p(i) %pb(i)=p(i-1)+krj(i)*(pb-p0); Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl)); Vw(i)=V-Vr(i); rr(i)=((3*Vr(i))/(4*pi))^(1/3); u(i)=((rr(i)-rr(i-1))/tj); mdw(i)=md0*((p(i)/p0)^(1/kl)); %Tw(i)=T0*((p(i)/p0)^((kl-1)/kl)); %Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); %mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i)); Vbw(i)=(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl)); Vbr(i)=V-Vbw(i); mdbw(i)=mdw(i-1)*(Vw(i-1)-Vrj(i))/Vbw(i); mdbr(i)=(mdr(i-1)*Vr(i-1)+mdw(i-1)*Vrj(i))/Vbr(i); Sl=150;Sr=150; Vxr(i)=((Sr*tj)/mdbr(i))*((2*p(i-1)*mdr(i-1)*(((p0/p(i-1))^(2/kr))-((p0/p(i-1))^((kr+1)/kr))))^(1/2)); Vxl(i)=((Sl*tj)/mdbw(i))*((2*p(i-1)*mdw(i-1)*(((p0/p(i-1))^(2/kl))-((p0/p(i-1))^((kl+1)/kl))))^(1/2)); pr(i)=pb(i)*(((Vbw(i)-(Vxl(i)))/Vbw(i))^(kl)); pl(i)=pb(i)*(((Vbr(i)-(Vxr(i)))/Vbr(i))^(kr)); p(i)=((Vbw(i)*((pr(i))^(1/kl))+((pl(i)/((pr(i)+pl(i))/2))^((kl-kr)/(kl*kr)))*Vbr(i)*((pl(i))^(1/kl)))/V)^kl; end end plot(t,p) |
回复此楼» 猜你喜欢
咱们一起用铁证分析2026国家社科基金中标与否
已经有21人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有5人回复
-
国基金的申报应该改成非等额制，评价高的钱多评价低的钱少，但是增加资助率
已经有8人回复
-
各位道友，我要去昆明玩几天，回来见。
已经有7人回复
-
欢迎发来filecode的Mz6后的代码验证其规律
已经有48人回复
-
filecode
已经有10人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有11人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有8人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有37人回复
-
好奇怪的filecode
已经有7人回复
pdl9527
专家顾问 (小有名气)
-
专家经验: +8 - 应助: 100 (初中生)
- 金币: 2110.3
- 红花: 23
- 帖子: 282
- 在线: 112小时
- 虫号: 1227333
- 注册: 2011-03-09
- 性别: GG
- 专业: 分离过程
- 管辖: 计算模拟
2楼2016-12-16 18:13:48
|
谢谢回复，我现在知道当p>pjc时是第318步，将选择条件改为while i<=318时进行选择计算（具体程序如下），但是现实结果出现In an assignment A(I) = B, the number of elements in B and I must be the same. Error in ==> jingque1 at 41 p(i)=p(i-1)+krj(i)*(pb-p0);%压力 实在想不出来这一步究竟错在哪里啊。 clear all p0=101.3;kp=7.9;T0=300;md0=1.3;u0=0.5;kz=4;%初始条件(压力、压缩系数、温度、密度、燃烧速度、湍流因子) Tbv=2221;Tbp=1867;%燃烧产物温度 V=5000;kl=1.4;kr=1.25;%屋室体积、冷气体与热气体比热容比 pb=p0*kp;%爆后压力 pjc=p0+5; t(1)=0;tj=0.001;%初始时间与时间间隔 Vr(1)=0;%初始燃烧产物体积 Vrj(1)=0;%单位时间内燃烧的气体体积 p(1)=p0;%初始压力 u(1)=u0;%初始速度 krj(1)=0;%初始燃尽率 rr(1)=0;%燃烧产物半径 mdw(1)=md0;%未燃气体密度 Tw(1)=T0;%未燃气体温度 i=1; %开始循环计算 while p(i)<=pb & p(i)>=p0 while i<318 i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25%正常传播速度--------------需要改进 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1);%单位时间内燃烧的体积 krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl));%燃尽率 p(i)=p(i-1)+krj(i)*(pb-p0);%压力 Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%燃烧产物体积 Vw(i)=V-Vr(i);%未燃产物体积 rr(i)=((3*Vr(i))/(4*pi))^(1/3);%燃烧产物半径 u(i)=((rr(i)-rr(i-1))/tj);%可见传播速度 mdw(i)=md0*((p(i)/p0)^(1/kl));%未燃气体密度 Tw(i)=T0*((p(i)/p0)^((kl-1)/kl));%未燃气体温度 Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i));%燃烧产物密度 end while i>=318 i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25%正常传播速度--------------需要改进 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1);%单位时间内燃烧的体积 krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl));%燃尽率 p(i)=p(i-1)+krj(i)*(pb-p0);%压力 pb(i)=p(i) %pb(i)=p(i-1)+krj(i)*(pb-p0);%压力 Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%燃烧产物体积 Vw(i)=V-Vr(i);%未燃产物体积 rr(i)=((3*Vr(i))/(4*pi))^(1/3);%燃烧产物半径 u(i)=((rr(i)-rr(i-1))/tj);%可见传播速度 mdw(i)=md0*((p(i)/p0)^(1/kl));%未燃气体密度 %Tw(i)=T0*((p(i)/p0)^((kl-1)/kl));%未燃气体温度 %Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); %mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i));%燃烧产物密度 Vbw(i)=(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%爆室内未燃气体体积 Vbr(i)=V-Vbw(i);%爆室内已燃气体体积 mdbw(i)=mdw(i-1)*(Vw(i-1)-Vrj(i))/Vbw(i);%爆室内未反应气体密度 mdbr(i)=(mdr(i-1)*Vr(i-1)+mdw(i-1)*Vrj(i))/Vbr(i);%爆室内燃烧产物密度 Sl=150;Sr=150; Vxr(i)=((Sr*tj)/mdbr(i))*((2*p(i-1)*mdr(i-1)*(((p0/p(i-1))^(2/kr))-((p0/p(i-1))^((kr+1)/kr))))^(1/2));%泄露的燃烧产物体积 Vxl(i)=((Sl*tj)/mdbw(i))*((2*p(i-1)*mdw(i-1)*(((p0/p(i-1))^(2/kl))-((p0/p(i-1))^((kl+1)/kl))))^(1/2));%泄露的未燃气体体积 pr(i)=pb(i)*(((Vbw(i)-(Vxl(i)))/Vbw(i))^(kl));%热气体过程压力 pl(i)=pb(i)*(((Vbr(i)-(Vxr(i)))/Vbr(i))^(kr));%冷气体过程压力 p(i)=((Vbw(i)*((pr(i))^(1/kl))+((pl(i)/((pr(i)+pl(i))/2))^((kl-kr)/(kl*kr)))*Vbr(i)*((pl(i))^(1/kl)))/V)^kl; end end plot(t,p) 请问这是为什么？ |
3楼2016-12-16 19:48:56
|
如果改成if语句，则可以计算但是还是只能算p<pjc之前的过程，请问是为什么，想了好久都不明白。 p0=101.3;kp=7.9;T0=300;md0=1.3;u0=0.5;kz=4;%初始条件(压力、压缩系数、温度、密度、燃烧速度、湍流因子) Tbv=2221;Tbp=1867;%燃烧产物温度 V=5000;kl=1.4;kr=1.25;%屋室体积、冷气体与热气体比热容比 pb=p0*kp;%爆后压力 pjc=p0+5; t(1)=0;tj=0.001;%初始时间与时间间隔 Vr(1)=0;%初始燃烧产物体积 Vrj(1)=0;%单位时间内燃烧的气体体积 p(1)=p0;%初始压力 u(1)=u0;%初始速度 krj(1)=0;%初始燃尽率 rr(1)=0;%燃烧产物半径 mdw(1)=md0;%未燃气体密度 Tw(1)=T0;%未燃气体温度 i=1; %开始循环计算 while p(i)<=pb & p(i)>=p0 if i<318 i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25%正常传播速度--------------需要改进 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1);%单位时间内燃烧的体积 krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl));%燃尽率 p(i)=p(i-1)+krj(i)*(pb-p0);%压力 Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%燃烧产物体积 Vw(i)=V-Vr(i);%未燃产物体积 rr(i)=((3*Vr(i))/(4*pi))^(1/3);%燃烧产物半径 u(i)=((rr(i)-rr(i-1))/tj);%可见传播速度 mdw(i)=md0*((p(i)/p0)^(1/kl));%未燃气体密度 Tw(i)=T0*((p(i)/p0)^((kl-1)/kl));%未燃气体温度 Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i));%燃烧产物密度 else i=i+1; t(i)=t(i-1)+tj; uz(i)=u0*(p(i-1)/p0)^0.25%正常传播速度--------------需要改进 Vrj(i)=(4*pi/3)*(((3*Vr(i-1)/(4*pi))^(1/3)+(kz*uz(i)*tj))^3)-Vr(i-1);%单位时间内燃烧的体积 krj(i)=(Vrj(i)/V)*((p(i-1)/p0)^(1/kl));%燃尽率 p(i)=p(i-1)+krj(i)*(pb-p0);%压力 pb(i)=p(i) %pb(i)=p(i-1)+krj(i)*(pb-p0);%压力 Vr(i)=V-(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%燃烧产物体积 Vw(i)=V-Vr(i);%未燃产物体积 rr(i)=((3*Vr(i))/(4*pi))^(1/3);%燃烧产物半径 u(i)=((rr(i)-rr(i-1))/tj);%可见传播速度 mdw(i)=md0*((p(i)/p0)^(1/kl));%未燃气体密度 %Tw(i)=T0*((p(i)/p0)^((kl-1)/kl));%未燃气体温度 %Tr(i)=Tbp+(Tbv-Tbp)*((p(i)-p0)/(pb-p0)); %mdr(i)=(p(i)*Tbv*md0)/(pb*Tr(i));%燃烧产物密度 Vbw(i)=(V-Vr(i-1)-Vrj(i))*((p(i-1)/p(i))^(1/kl));%爆室内未燃气体体积 Vbr(i)=V-Vbw(i);%爆室内已燃气体体积 mdbw(i)=mdw(i-1)*(Vw(i-1)-Vrj(i))/Vbw(i);%爆室内未反应气体密度 mdbr(i)=(mdr(i-1)*Vr(i-1)+mdw(i-1)*Vrj(i))/Vbr(i);%爆室内燃烧产物密度 Sl=150;Sr=150; Vxr(i)=((Sr*tj)/mdbr(i))*((2*p(i-1)*mdr(i-1)*(((p0/p(i-1))^(2/kr))-((p0/p(i-1))^((kr+1)/kr))))^(1/2));%泄露的燃烧产物体积 Vxl(i)=((Sl*tj)/mdbw(i))*((2*p(i-1)*mdw(i-1)*(((p0/p(i-1))^(2/kl))-((p0/p(i-1))^((kl+1)/kl))))^(1/2));%泄露的未燃气体体积 pr(i)=pb(i)*(((Vbw(i)-(Vxl(i)))/Vbw(i))^(kl));%热气体过程压力 pl(i)=pb(i)*(((Vbr(i)-(Vxr(i)))/Vbr(i))^(kr));%冷气体过程压力 p(i)=((Vbw(i)*((pr(i))^(1/kl))+((pl(i)/((pr(i)+pl(i))/2))^((kl-kr)/(kl*kr)))*Vbr(i)*((pl(i))^(1/kl)))/V)^kl; end end plot(t,p) |
4楼2016-12-16 20:03:25
pdl9527
专家顾问 (小有名气)
-
专家经验: +8 - 应助: 100 (初中生)
- 金币: 2110.3
- 红花: 23
- 帖子: 282
- 在线: 112小时
- 虫号: 1227333
- 注册: 2011-03-09
- 性别: GG
- 专业: 分离过程
- 管辖: 计算模拟
5楼2016-12-16 21:21:38
pdl9527
专家顾问 (小有名气)
-
专家经验: +8 - 应助: 100 (初中生)
- 金币: 2110.3
- 红花: 23
- 帖子: 282
- 在线: 112小时
- 虫号: 1227333
- 注册: 2011-03-09
- 性别: GG
- 专业: 分离过程
- 管辖: 计算模拟
6楼2016-12-16 21:29:35
mygt_hit
专家顾问 (职业作家)
-
专家经验: +362 - 应助: 438 (硕士)
- 贵宾: 0.019
- 金币: 19579.9
- 散金: 5130
- 红花: 135
- 沙发: 3
- 帖子: 4939
- 在线: 990.5小时
- 虫号: 1489764
- 注册: 2011-11-13
- 性别: GG
- 专业: 结构工程
- 管辖: 土木建筑
7楼2016-12-17 13:01:04
