>	restart:interface(warnlevel=0):#"m01_p11"

Dentro de un cilindro vertical de 30 litros y 1 dm2 de seccin, sellado por ambos extremos, existe un émbolo de 4 kg que separa dos masas iguales de aire. Inicialmente la presin en la parte superior es la atmásférica, y entonces se suelta desde la parte superior del cilindro una pesa de 8 kg. Se pide:

a) Presin inicial en la parte inferior, altura del émbolo, tiempo que tarda en chocar la pesa y velocidad con la que llega y con la que empieza a moverse el émbolo.

b)Altura final del émbolo, presiones finales y flujo de calor hasta el equilibrio final.

c)Altura del émbolo y presiones tras el equilibrio mecánico suponiendo que toda la disipacin se concentra en el émbolo.

d)Altura mínima alcanzada por el émbolo y presiones en ese instante.

Datos:

>	read"../therm_eq.m":read"../therm_proc.m":with(therm_proc):

>	su:="Aire":dat:=[Vt=0.03*m_^3,A=0.01*m_^2,mE=4*kg_,mP=8*kg_,p11=subs(Const,p0)];

Esquema:

>

>

Ecs. const.:

>	eqET:=subs(eq1_11,eq1_12);eqEE:=eq1_16;gdat:=get_gas_data(su):dat:=[op(dat),Const,gdat,SI2,SI1]:

a) Presin inicial en la parte inferior, altura del émbolo, tiempo que tarda en chocar la pesa y velocidad con la que llega y con la que empieza a moverse el émbolo.

>

ecBF:=m[EP]*diff(z(t),t,t)=(p[2]-p[1])*A-m[EP]*g-Ff;p21:=p11+mE*g/A;p21_:=subs(dat,SI1,p21):p21__:=evalf(%/(1000*Pa_/kPa));L_:=subs(dat,Vt/A);eqET1:=subs(V=Vt-z1*A,p=p11,T=T0,eqET);eqET1_:=subs(dat,eqET1):eqET2:=subs(V=z1*A,p=p21,T=T0,eqET);eqET2_:=subs(dat,eqET2):sol1:=subs(dat,solve({eqET1_,eqET2_},{m,z1})):dat:=[op(sol1),L=L_,op(dat)]:z1_:=subs(sol1,z1);m__:=subs(sol1,m);

>	eqMov:=L-z1=(1/2)*g*t^2;t01_:=evalf(solve(subs(dat,eqMov),t)[1]):t01=evalf(%,2);vel1:=g*t;vel1_:=evalf(subs(t=t01_,dat,vel1)):vel1=evalf(%,2);eqCh:=mP*vel1=(mP+mE)*vel11;vel11_:=evalf(solve(subs(t=t01_,dat,eqCh),vel11)):vel11=evalf(%,2);

b)Altura final del émbolo, presiones finales y flujo de calor hasta el equilibrio final.

>

p22:=solve(subs(V=A*z2,T=T0,eqET),p);p12:=solve(subs(V=Vt-A*z2,T=T0,eqET),p);eqBF:='p22=p12+(mP+mE)*g/A';z2_:=evalf(solve(subs(dat,SI2,SI2,eqBF),z2))[2];dat:=['z2'=z2_,op(dat)]:p22_:=subs(dat,SI2,SI1,p22):'p22'=evalf(%/(1e3*Pa_/kPa_));p12_:=subs(dat,SI2,SI1,p12):'p12'=evalf(%/(1e3*Pa_/kPa_));Q12:=mP*g*(L-z2)+mE*g*(z1-z2);Q12_:=subs(dat,SI1,Q12);

c)Altura del émbolo y presiones tras el equilibrio mecánico suponiendo que toda la disipacin se concentra en el émbolo.

>

p23:=solve(subs(V=A*z3,T=T23,eqET),p);p13:=solve(subs(V=Vt-A*z3,T=T13,eqET),p);T23:='T0*(z1/z3)^(gamma-1)';T13:='T0*((L-z1)/(L-z3))^(gamma-1)';eqBF:='p23=p13+(mP+mE)*g/A';z3_:=evalf(fsolve(subs(dat,SI0,eqBF),z3=0..3))*m_;dat:=[op(dat),'z3'=z3_]:p23_:=subs(dat,SI2,SI1,p23):'p23'=evalf(%/(1e3*Pa_/kPa_));p13_:=subs(dat,SI2,SI1,p13):'p13'=evalf(%/(1e3*Pa_/kPa_));T23_:=evalf(subs(dat,T23));T13_:=evalf(subs(dat,T13));

d)Altura mínima alcanzada por el émbolo y presiones en ese instante.

>

p24:=solve(subs(V=A*z4,T=T24,eqET),p);p14:=solve(subs(V=Vt-A*z4,T=T14,eqET),p);T24:='T0*(z1/z4)^(gamma-1)';T14:='T0*((L-z1)/(L-z4))^(gamma-1)';eqBE4:=(mP+mE)*vel11^2/2+(mP+mE)*g*(z1-z4)=m*c[v]*(T24-T0+T14-T0);z4_:=evalf(fsolve(subs(vel11=vel11_,dat,SI0,eqBE4),z4=0..3)*m_):z4=evalf(%,3);dat:=['z4'=z4_,op(dat)]:p24_:=subs(dat,SI2,SI1,p24):'p24'=evalf(%/(1e3*Pa_/kPa_));p14_:=subs(dat,SI2,SI1,p14):'p14'=evalf(%/(1e3*Pa_/kPa_));T24_:=evalf(subs(dat,T24));T14_:=evalf(subs(dat,T14));

e) Dinámica

>	coe_:=30:p1:=p11*((L-z1)/(L-z(t)))^gamma;p2:='p21'*(z1/z(t))^gamma;deq1:=diff(z(t),t)=zp(t);deq2:=(mP+mE)*diff(zp(t),t)=(p2-p1)*A-(mP+mE)*g-coe*zp(t);ic1:=z(0)=z1;ic2:=zp(0)=-vel11;deq1_:=deq1:deq2_:=subs(coe=coe_,dat,SI0,deq2);ic1_:=subs(dat,SI0,ic1):ic2_:=subs(vel11=vel11_,SI0,ic2):

>	dsol1:=dsolve({deq1_,deq2_,ic1_,ic2_},{z(t),zp(t)},numeric);

>	with(plots):odeplot(dsol1,[t,z(t)],0..2,numpoints=500,labels=["t(seg)","z(m)"],colour=black);odeplot(dsol1,[t,zp(t)],0..2,numpoints=500,labels=["t(seg)","v(m/s)"],colour=black);odeplot(dsol1,[z(t),zp(t)],0..1,numpoints=100,labels=["z(m)","v(m/s)"],colour=black);

Lineal

>

acel:=convert(series(subs(z(t)=z1+Dz,coe=0,rhs(deq2)),Dz=0,2),polynom);acelmod:=(mP+mE)*(-(2*Pi/T)^2)*(Dz-Dz0);eqeq:='acel-acelmod';eqeq:=collect(acel-acelmod,Dz)=0:T_:=subs(dat,SI0,solve(coeff(lhs(eqeq),Dz,1),T)[1])*s_:Dz0_:=solve(coeff(lhs(eqeq),Dz,0),Dz0):'T'=evalf(T_,2);Dz0=evalf(subs(T=T_,dat,Dz0_),2);

>