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);