function wpy531
clear all
clc

[t,U]=ode45(@fly1,[0:0.0001:0.005],0.79);
figure(1)
plot(t,U)
title('Droplet Velocity VS Time')
xlabel('Time/s')
ylabel('Droplet Velocity/m/s')

[t,T]=ode45(@temperaturel,[0:0.0001:0.005],1173,[],U(end));
figure(2)
plot(t,T)
title('Droplet Velocity VS Time')
xlabel('Time/s')
ylabel('Droplet Velocity/m/s')

function y=fly1(t,U)
a=1.1616;                  %µªÆøÃܶÈ
b=1.873e-5;                %µªÆøð¤¶È
d=1071e-6;                    %ÈÛµÎÖ±¾¶
c=2385;                    %ÈÛµÎÃܶÈ
g=9.8;                     %ÖØÁ¦¼ÓËÙ¶È
f=d*a*U/b;                 %ReÊý
h=0.28+6/f.^0.5+21/f;      %ÍÏҷϵÊý
y=((1-a/c)*g-3*(0.28+6/f.^0.5+21/f)*a*U.^2/(4*d*c))/(1+0.5*a/c);

function x=temperaturel(t,T,U)% ÒºÏàÀäÈ´
rhod=2368;                % rho-d, ÈÛµÎÃܶȣ¬µ¥Î»Îªkg*m-3,-0.28*(T-933)
cpd=1086;                %cpd, È۵αÈÈÈÈÝ£¬µ¥Î»ÎªJm-3K-1
zhj=1.102e-3;            %dd,ÈÛµÎÖ±¾¶£¬µ¥Î»m
d=0.6;                  %È۵κڶȣ¬
f=5.6032e-8;             %Stefan-Boltzmann³£Êý£¬µ¥Î»W*m-3*K-4
g=293;                  %T0£¬»·¾³Î¶ȣ¬µ¥Î»K
ug=2.217e-5;              %ë²ÆøÕ³¶È
cpg=520.3;                 %ë²ÆøÈÈÈÝ
kg=0.01712;                %kg,µªÆøÈȵ¼ÂÊ£¬µ¥Î»W*m-1*K-1
rhog=1.784;               %rhog, ÆøÌåÃܶȣ¬µ¥Î»kg*m-3
sudu=3.2+7.27*t;                   %V£¬ÈÛµÎËٶȣ¬µ¥Î»m*s-1
h=kg/zhj*(2.0+0.6*(rhog*zhj*U/ug)^0.5*(ug*cpg/kg)^(1/3));  %h, È۵ζÔÁ÷´«ÈÈϵÊý£¬W/(m2*K)
x=-6*(0.6*(5.6032e-8)*(T^4-7.3701e+009)+h*(T-293))/(rhod*cpd*zhj);