Help me to solve the function

2020 9 29
1 个回答
回答已采纳
更新时间：2020 10 19
1 次查看（30 天）

 采纳的回答

KSSV
KSSV 2020-9-29

0 个投票

clear; clc;
wc=0.4;
pc=3.16;
T=293;
rwc3s=0.586;
rwc3a=0.172;
mc=0.3;
vc=1/((wc*pc)+1);
ro=(3*vc/(4*pi))^1/3;
yg=0.25;
yw=0.15;
RH=0.88;
b1=1231;
b2=7579;
B293=0.3*10^-10;
C=5*10^7;
ER=5364;
De293=((rwc3s*100)^2.024)*(3.2*10^-14);
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);
v=2.15;
cw=((RH-0.75)/0.25)^3;
pw=1;
B=B293*exp(-b1*(1/T-1/293));
kr=kr293*exp(-ER*(1/T-1/293));
da=360;
dt=1;
Nn=(da)*(1/dt);
timestep=0;
time=1:Nn;
time2=1:Nn+1;
% r_initial=16;
r_initial=[ 0.5 2 8 16 32];
m = length(r_initial) ;
n1 = length(time) ; n2 = length(time2) ;
Alfa = zeros(m,n1) ;
RT_INC1 = zeros(m,n1) ;
RT_INC2 = zeros(m,n1) ;
for i = 1:length(r_initial)
ro_init=r_initial(i)*0.001;
L=((4*pi*(wc*pc/pw+1)/3)^(1/3))*ro_init;
for BB=1:Nn
if (BB==1)
rt_inc(1)=ro_init(1)-0.00001; %boundary condition for unhydrated
else
rt_inc(1)=ro_init(1)-0.00001;
end
rt_inc(Nn+1)=0;
if (BB==1)
Rt_inc(1)=ro_init(1)+0.00001; %boundary condition for hydrated
else
Rt_inc(1)=ro_init(1)+0.00001;
end
Rt_inc(Nn+1)=0;
end
for M=1:Nn
if Rt_inc(M)/L<=ro
Sr(M)=0;
elseif (Rt_inc(M)/L>=ro)&&(Rt_inc(M)/L<0.5)
Sr(M)=4*pi*(Rt_inc(M)/L)^2;
elseif (0.5<=Rt_inc(M)/L)&&(Rt_inc(M)/L<0.5*(2^0.5))
Sr(M)=(4*pi*(Rt_inc(M)/L)^2)-(6*pi*(1-(0.5/(Rt_inc(M)/L))));
elseif (Rt_inc(M)/L>=0.5*(2^0.5)) && (Rt_inc(M)/L<0.5*(3^0.5))
fun = @(y,x) 8*(Rt_inc(M)/L)./(sqrt((Rt_inc(M)/L)^2-(x.^2)-(y.^2)));
ymin=sqrt((Rt_inc(M)/L)^2-0.5);
xmin=@(x) sqrt((Rt_inc(M)/L)^2-0.25-x.^2);
Sr(M)=integral2(fun,ymin,0.5,xmin,0.5);
else
Sr(M)=0;
end
cst(M)=Sr(M)/(4*pi*Rt_inc(M)^2);
alfa(M)=1-(rt_inc(M)/ro_init)^3;
kd(M)=(B/(alfa(M)^1.5))+C*(Rt_inc(M)-ro_init)^4;
De(M)=De293*(log(1/alfa(M)))^1.5;
rt_inc(M+1)=rt_inc(M)-(dt*((pw*cst(M)*cw)/((yw+yg)*pc*rt_inc(M)^2))*1/((1/(kd(M)*rt_inc(M)^2))+(((1/rt_inc(M))-(1/Rt_inc(M)))/De(M))+(1/(kr*rt_inc(M)^2))));
Rt_inc(M+1)=((v-1)*((rt_inc(M)-rt_inc(M+1))/dt))*dt+Rt_inc(M);
end
Alfa(i,:) = alfa ;
RT_INC1 = Rt_inc ;
RT_INC2 = rt_inc ;
end
figure(1)
plot(time,Alfa)
legend
xlabel('Time (Hours)')
xlim([0 da])
grid on
figure (2)
plot(time2,RT_INC1,time2,RT_INC2),grid;
legend
yline(ro_init,'-','Initial Radius Particles');
legend('R_t','r_t','r_o')
%title('xx')
xlabel('Time (Hours)')
ylabel('Radius of Particles (mm)')
xlim([0 da])
grid on

3 个评论
显示 1更早的评论 隐藏 1更早的评论

Kevin Isakayoga
Kevin Isakayoga 2020-9-29
Thank you very much !
Kevin Isakayoga
Kevin Isakayoga 2020-9-29
By the way, I think the last graph is not the same like I want, because its only depends on the final value which in this case is 32.
r_initial=[ 0.5 2 8 16 32];
Could you help me to make each of the r_initial value influence to become a graph? I already tried like this,
RT_INC1(i,:) = Rt_inc ;
RT_INC2(i,:) = rt_inc ;
the value become messy.
figure (2)
plot(time2,RT_INC1,time2,RT_INC2),grid;
yline(ro_init,'-','Initial Radius Particles');
legend('R_t','r_t','r_o')
xlabel('Time (Hours)')
ylabel('Radius of Particles (mm)')
xlim([0 da])
grid on
Thank you very much in advance! Looking forward to your reply.
KSSV
KSSV 2020-9-29
It is changing.....you need to check your formula..

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Frequency-Domain Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by