气动仿真作业2008-5-14

1、气动仿真作业2008-5-14张业明 BY07031301 题目描述条件：给定一个体积为30L的容器，对其充气，声速流导C=0.2dm3/(s·bar)；有效截面积Se=1mm3；临界压力比b=0.4；入口压力Ps=0.5MPa；室温为20大气压力为1bar。问题：（1） 多变过程（n=1，n=1.4），求P（t）和（t）；（2） 传热率h=20W/m2·K。2 问题分析（1）多变过程多变过程中满足的方程式如下： (1)这里，n是多变指数。等温过程时，n=1；绝热过程时，n=k=1.4。对（1）进行微分，可得(2)质量流量G的公式: (3)质量流量G的单位为kg/s。式中：

2、标准状态下的体积流量dm3/s(ANR)C ：声速流导dm3/(s·bar)P1：上流绝对压力barP2：下流绝对压力bar ：上流空气的绝对温度Kb ：临界压力比：标准状态下的密度为1.185kg/m3声速流导C和有效截面积Se尽管单位不同，但都是表示流通能力的参数。两者之间存在如下的换算关系：Semm2=5.03×Cdm3/(s·bar) 5×Cdm3/(s·bar)(4)由方程（4）可以求得声速流导C,再由（3）可以得到质量流量G，将G代入（2）得到关于t的微分方程。求解得到。（2）传热过程（热传率一定）根据能量守恒定律，容腔内空气的内能

3、的变化量dU取决于流进或流出的空气的焓、外界的热交换量（5）容腔壁通常由金属制成，金属的热容量比空气的热容量大很多，金属内壁的温度可假设为不受空气充放气影响，与室温相同。这样，容腔内空气与容腔内壁的热交换可用下式表示（6）这里，是空气与壁面的热传率，是容腔内壁的表面积。将（6）代入（5），可得（7）这里，是上流的温度，充气的时候是气源温度，放气的时候是容腔内的空气温度。将（7）左侧展开后整理得到（8）由压力微分方程（状态方程） 对时间微分 (9)由（8）和（9）联立后整理，得到如下的微分方程组（10）是上流的温度，充气的时候是气源温度，放气的时候是容腔内的空气温度；是室温。将G代入（10）求得

4、。3 问题的Matlab仿真求解（1）多变过程的仿真程序及仿真曲线主求解程序solve1SYS.mfunction solve1SYS()p=1;t0=293;options = odeset('RelTol',1e-4,'AbsTol',1e-4 1e-4);T,Y = ode45(solve1,0 5,p t0,options);subplot(2,1,1);plot(T,Y(:,1),'r-');xlabel('t/s'),ylabel('P/bar');title('P(t)');legen

5、d('P(t)');subplot(2,1,2);plot(T,Y(:,2),'b-');xlabel('t/s'),ylabel('/K');title('(t)');legend('(t)');问题描述文件solve1.mfunction dy=solve1(t,y)%空气常数：R=287J/(kg.K)%容腔体积：v=30L%声速流导：C=0.2dm3/(s.bar)%有效截面积：Se=5.03*C=1.006 k=1.4; %绝热%k=1; %等温C=0.2; %dm3/(s.bar)R=2

6、87; %J/(kg.K)V=30; %Lp1=6; %bart1=293; %Kb=0.4; %density=1.185e-3 ; %kg/dm*3 dy=zeros(2,1);if(y(1)/p1<=b) G=density*C*p1*sqrt(293/t1);elseif(y(1)/p1<=1)G=density*C*p1*sqrt(293/t1)*sqrt(1-(y(1)/p1-b)/(1-b)2);else G=0;enddy(1)=k*G*R*y(2)/V;dy(2)=(k-1)*G*R*y(2)2/y(1)/V; 绝热条件下的仿真结果：等温条件下的仿真结果：(2)传

7、热过程的仿真程序及仿真曲线主求解程序solve2SYS.mfunction solve2SYS()p=1;t0=293;options = odeset('RelTol',1e-4,'AbsTol',1e-4 1e-4);T,Y = ode23t(solve2,0 5,p t0,options);subplot(2,1,1);plot(T,Y(:,1),'r-');xlabel('t/s'),ylabel('P/bar');title('P(t)');legend('P(t)');s

8、ubplot(2,1,2);plot(T,Y(:,2),'b-');xlabel('t/s'),ylabel('/K');title('(t)');legend('(t)');问题描述文件solve2.mfunction dy=solve2(t,y)%空气常数：R=287J/(kg.K)%容腔体积：v=30L%声速流导：C=0.2dm3/(s.bar)%有效截面积：Se=5.03*C=1.006 K=1.4; %绝热%k=1; %等温C=0.2; %dm3/(s.bar)R=287; %J/(kg.K)V=30;

9、%Lp1=6; %bart1=293; %Kb=0.4; %density=1.185e-3 ; %kg/dm*3Cv=R/(K-1);Cp=Cv+R;h=20e3;Sh=4*pi*(V/(4*pi/3)(2/3);ta=293; dy=zeros(2,1);if(y(1)/p1<=b) G=density*C*p1*sqrt(293/t1);elsey(1)/p1<=1)G=density*C*p1*sqrt(293/t1)*sqrt(1-(y(1)/p1-b)/(1-b)2);elseG=0;enddy(1)=R/Cv/V*(G*Cp*t1+h*Sh*(ta-y(2);dy(2

10、)=R*y(2)/Cv/y(1)/V*(G*Cp*t1-G*Cv*y(2)+h*Sh*(ta-y(2); 热传率一定的条件下的仿真结果：4 容腔充放气的无因次数学模型4.1 基准量(1)基准压力 充气时定为供气压力；放气时定为容腔的初始压力。(2)基准温度 定为外界环境温度，即室温。(3)基准流量Gmax 定为基准压力下的声速流的流量，也是充气或放气过程中的最大流量(11)这里，=1.185kg/m3,=293K,C是充气或放气通道的声速流导。（4）基准时间充气时的基准时间是以基准流量将容腔从完全真空充到基准压力所需的时间；放气时的基准时间是以基准流量将容腔从基准压力放气到完全真空所需的时间。

11、其计算表达式如下(12)4.2 无因次数学模型4.2.1 容腔充气首先，对容腔内压力、温度和时间进行无因次化(13)(14)(15)通常，供给空气的温度就是室温，这里就假设两者相同。这样，充入容腔的空气流量可表示为 (16)式（10）的压力微分式可写成（17）将其无因次化，可得(18)其中，热平衡时间常数 （19）其物理意义为通过传热将容腔内全部空气具有的内能传递给绝对零度环境所需的时间。（20）将（19）、（20）代入（18）得（21）同样对（10）的温度微分式进行无因次化，可得（22）（23）将（21）和（22）联立，构成压力和温度的无因次微分方程组，即容腔充气的无因次数学模型（24）4.

12、2.2容腔放气从容腔放气的时候，流出容腔的空气流量可表示为（25）注意，在放气时，即使是声速流，因为空气流量与容腔内温度有关，所以与充气不同，不是定值。放气时上流温度为容腔内温度，式（10）可以写成（26）与充气同样将其无因次化，可得（27）4.3 容腔充放气的无因次仿真程序和仿真结果充气仿真主程序文件solve3SYS.mfunction solve3SYS() Ps=6; % bar P0=1; %环境压力 ndm_P0=P0/Ps; Q1=293; % ANR Qa=293; % 室温 ndm_Q0=1; options = odeset('RelTol',1e-4,&#

13、39;AbsTol',1e-4 1e-4);T,Y = ode23t(solve3,0 5,ndm_P0 ndm_Q0,options); subplot(3,1,1);hold onplot(T,Y(:,1),'r-'); grid on;xlabel('t/s')title('P*(t)');legend('P*(t)'); subplot(3,1,2);hold onplot(T,Y(:,2),'b-');grid on;xlabel('t/s')title('*(t)'

14、;);legend('*(t)'); b=0.4;m,n=size(Y(:,1);G=zeros(size(Y(:,1);for i=1:m if(Y(i,1)<=b) G(i,1)=1; elseif(Y(i,1)<1) G(i,1)=sqrt(1-(Y(i,1)-b)/(1-b)2); else G(i,1)=0; end end subplot(3,1,3);hold onplot(T,G,'g-');grid on;xlabel('t/s')title('G*(t)');legend('G*(t)

15、9;);问题描述文件solve3.mfunction dy=solve3(t,y)%空气常数：R=287J/(kg.K)%容腔体积：v=30L%声速流导：C=0.2dm3/(s.bar)%有效截面积：Se=5.03*C=1.006 b=0.4; V=30; %单位dm3 C=0.2; %声速流导 p0=1.185e-003; %密度 kg/dm3 注意单位 Ps=6; % bar P0=1; %环境压力 ndm_P0=P0/Ps; h=20; Sh=4*pi*(V/(4*pi/3)(2/3) ; % 46.6908 R=287; Q1=293; % ANR Qa=293; % 室温 ndm_Q

16、0=1; y0=ndm_P0 ndm_Q0; K=1.4; Cv=R/(K-1); Cp=Cv+R; Tp=V/(R*Qa*P0*C)*sqrt(Qa/Q1); Th=Cv*Ps*V/(h*Sh*R*Qa); % Th=m*Cv/(h*Sh) 注意单位 %Ka=Tp/Th % Ka=Tp/Th=1.0846 %Ka=2 Ka=10; dy=zeros(2,1); if(y(1)<=b) G=1;elseif(y(1)<1) G=sqrt(1-(y(1)-b)/(1-b)2);else G=0;end dy(1)=K*G+Ka*(1-y(2);dy(2)=y(2)/y(1)*(K-y

17、(2)*G+Ka*(1-y(2); 仿真结果：无因次充气ka=1.0846.jpg无因次充气ka=1.0846 ka=2.jpg无因次充气ka=1.0846 ka=2 ka=10.jpg放气仿真主程序solve4SYS.mfunction solve4SYS() Ps=6; % bar P0=1; %环境压力 ndm_P0=P0/Ps; Q1=293; % ANR Qa=293; % 室温 ndm_Q0=1; options = odeset('RelTol',1e-4,'AbsTol',1e-4 1e-4);T,Y = ode23t(solve4,0 5,1

18、ndm_Q0,options); subplot(3,1,1);hold onplot(T,Y(:,1),'r-'); grid on;xlabel('t/s')title('P*(t)');legend('P*(t)'); subplot(3,1,2);hold onplot(T,Y(:,2),'b-');grid on;xlabel('t/s')title('*(t)');legend('*(t)'); b=0.4;m,n=size(Y(:,1);G=zeros(

19、size(Y(:,1);for i=1:m if(ndm_P0/Y(i,1)<=b) G(i,1)=Y(i,1)/sqrt(Y(i,2); elseif(ndm_P0/Y(i,1)<1) G(i,1)=Y(i,1)/sqrt(Y(i,2)*sqrt(1-(ndm_P0/Y(i,1)-b)/(1-b)2); else G(i,1)=0; end end subplot(3,1,3);hold onplot(T,G,'g-');grid on;xlabel('t/s')title('G*(t)');legend('G*(t)');问题描述文件solve4.mfunction dy=solve4(t,y)%空气常数：R=287J/(kg.K)%容腔体积：v=30L%声速流导：C=0.2dm3/(s.bar)%有效截面积