why I do not receive smooth transition and smooth lines

1 view (last 30 days)
hello dear friends.
could you please explain me why i do not receive smooth transition like below graph
the graph is taken from
clear
close all, clc
ca = 1014 ; % Specific Heat of Air, 𝑐a, J/kg K
ci = 2050 ; % Specific Heat of Ice, 𝑐𝑖, J/kg K
cw = 4218 ; % Specific Heat of Water, 𝑐𝑤, J/kg K
e0 = 27.03 ; % Sublimation and Evaporation Constant, 𝜀0 Pa/K
LWC = 0.001; % liquid water content, or G, kg/m3
Haw = 500 ; % Heat Transfer Coefficient, between air and water, W/m2 K
His = 999 ; % Heat Transfer Coefficient, between ice and substrate, W/m2 K
ki = 2.18 ; % thermal conductivity of ice in W/m K
kw = 0.571 ; % thermal conductivity of water in W/m K
Le = 13.4 ; % Lewis number, ka / (ca * mu_a)
LF = 3.344 * 1e5 ; % Latent Heat of Fusion, 𝐿𝐹, J/kg
LE = 2.26 * 1e6 ; % Latent Heat of Evaporation, 𝐿𝐸, J/kg
p0 = 6.3*1e4 ; % ambient pressure, in Pa
v = 90 ; % Freestream velocity
beta = 0.55 ; % Collection Efficiency, β
rho_g = 917; % Glaze Ice density, kg/m3
rho_r = 880; % Rime Ice density, kg/m3
rho_w = 1000; % Water density, kg/m3
X = 11 ; % Evaporation Coefficient
r = 0.9 ; % Recovery constant, used in Qa
lamda = 1.4 ; % Heat Capacity of Air
R = 287.05 ; % Gas Constant, 𝑅, J/kg K
Tf = 273.15 ; % freezing temperature in kelvin
cp = 1012 ; % specific heat of air J/(g.k), used in Qa
step=0.1; % step size
time=0:step:120; % 1201 steps
disp('*********** when Ta = 263.15 K *********** ')
*********** when Ta = 263.15 K ***********
Ta = 263.15;
Ts=Ta;
% Energy terms
Qa=r*Haw*((v)^2/(2*cp));
Qk=LWC*beta*((v^3)/2);
dBdt=((LWC * beta * v )/rho_r );
QL=rho_r*LF*dBdt;
denominator_Bg = LWC*beta*v*LF + (Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta);
Bg = (ki*(Tf-Ts)) / denominator_Bg; % Glaze ice thickness
Bgg=Bg*1e3 % to convert to mm
Bgg = 2.5668
tg = (rho_r /(LWC*beta*v))*Bg ; % time when first glaze appears
tg=str2num(num2str(tg,'%1.1f'))
tg = 45.6000
disp([ 'The glace ice first appears at ', num2str(tg,'%1.1f'),'s , and the thickness of the glaze ice is = ',num2str(Bg*1000),'mm'])
The glace ice first appears at 45.6s , and the thickness of the glaze ice is = 2.5668mm
% after finding the transition time, rime ice maximum thickness could be plotted
% find index of tg within 120 s
tg_index=find(abs(time-tg)<1e-3);
% use tg_index to create points between 0 and Bg, this is to plot the line up to transition.
% Br=linspace(0,(Bg*1000),tg_index);
tr=0:step:tg;
for i=1:length(tr)
Br(i)=((LWC*v*beta)/rho_r)*tr(i)*1e3;
end
% initial condition
t(1)=tg;
B(1)=Bg;
% Define the ODE function handle
% equation 4 from extended messinger
% substitute 32 and 45 equations after differentiate them into eqn 4
aa= 1/(rho_g*LF); % coefficent of dB/dt in eqn 4
bb= (ki * (Tf-Ts)/B); % dT/dz in equation 4
cc= ((kw*rho_w)*((Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta)))/((kw*rho_w)+(((Haw+v*beta)*(t-tg))-(rho_g*(B-Bg)))*(Haw+LWC*beta*v*cw+X*e0)); % dTheta/dz in eqn 4
% h=(((v+LWC*beta)*(t-tg))/rho_w)-((rho_g/rho_w)*(B-Bg));
% cc= (((Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta)))/(kw+h*(Haw+LWC*beta*v*cw+X*e0)); % dTheta/dz in eqn 4
f = @(t,B) aa*(bb-cc); % change the function as you desire
% RK4 loop
d=length(time)-length(tr); % start from transition time
for i = 1:d-1
% % % % % update x
t(i+1)= t(i)+step;
% % % % % update y
% k1 = step*(f(t(i) ,B(i) ));
% k2 = step*(f(t(i)+0.5*step ,B(i)+0.5*step*k1 ));
% k3 = step*(f(t(i)+0.5*step ,B(i)+0.5*step*k2 ));
% k4 = step*(F_xy(t(i)+ step ,B(i)+ step*k3 ));
% B(i+1) = B(i) + (step)*(k1+2*k2+2*k3);
k1 = step*f(t(i) , B(i)); k1 = double(k1);
k2 = step*f(t(i)+step/4 , B(i)+k1/4); k2 = double(k2);
k3 = step*f(t(i)+3/8*step , B(i)+3/32*k1+9/32*k2); k3 = double(k3);
gamma = [16/135, 0, 6656/12825];
K = [k1, k2, k3];
B(i+1) = B(i) + sum(K.*gamma);
end
t_263=[tr t];
B_263=[Br (B*1e3)];
% figure
% hold on
% grid on
% plot(t_263,B_263,'r');
% title(' Ice and water growth for Ta=263.15 K and Ta=270 K; Ts=Ta in both cases ');
% legend('B@ 263.15K','Position',[0.2, 0.78, .1, .1])
% xlabel(' Time(s) '); ylabel(' B,h (mm) ');
disp(' ')
disp('*********** when Ta = 270 K *********** ')
*********** when Ta = 270 K ***********
Ta = 270;
Ts=Ta;
% Energy terms
Qa=r*Haw*(v^2/(2*cp));
Qk=LWC*beta*((v^3)/2);
% dBdt=((LWC * beta * v )/rho_r );
% QL=rho_r*LF*dBdt;
denominator_Bg = LWC*beta*v*LF + (Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta);
Bg = (ki*(Tf-Ts)) / denominator_Bg; % Glaze ice thickness
Bgg=Bg*1e3 % to convert to mm
Bgg = 0.4463
tg = (rho_r /(LWC*beta*v))*Bg; % time when first glaze appears
tg=str2num(num2str(tg,'%1.1f'))
tg = 7.9000
disp([ 'The glace ice first appears at ', num2str(tg,'%1.1f'),'s , and the thickness of the glaze ice is = ',num2str(Bg*1000),'mm'])
The glace ice first appears at 7.9s , and the thickness of the glaze ice is = 0.44635mm
% after finding the transition time, rime ice maximum thickness could be plotted
tg_index=find(abs(time-tg)<1e-3);
Br=linspace(0,(Bg*1000),tg_index);
tr=0:step:tg;
% initial condition
tt(1)=tg;
BB(1)=Bg;
% Define the ODE function handle
% equation 4 from extended messinger
% substitute 32 and 45 equations after differentiate them into eqn 4
aa= 1/(rho_g*LF); % coefficent of dB/dt in eqn 4
bb= (ki * (Tf-Ts)/BB); % dT/dz in equation 4
cc= ((kw*rho_w)*((Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta)))/(kw*rho_w+((v+LWC*beta)*(tt-tg)-(rho_g*(BB-Bg)))*(Haw+LWC*beta*v*cw+X*e0)); % dTheta/dz in eqn 4
% h=(((v+LWC*beta)*(tt-tg))/rho_w)-((rho_g/rho_w)*(BB-Bg));
% cc= (((Qa+Qk)-(Haw+LWC*beta*v*cw+X*e0)*(Tf-Ta)))/(kw+h*(Haw+LWC*beta*v*cw+X*e0)); % dTheta/dz in eqn 4
F_xy = @(tt,BB) aa*(bb-cc); % change the function as you desire
% RK4 loop
d=length(time)-length(tr); % to start from transition time
for i = 1:d-1
% update x
tt(i+1)= tt(i)+step;
% update y
k1 = step*f(tt(i) , BB(i)); k1 = double(k1);
k2 = step*f(tt(i)+step/4 , BB(i)+k1/4); k2 = double(k2);
K = [k1, k2];
gamma = [16/135, 0.2];
BB(i+1) = BB(i) + sum(K.*gamma);
end
t_270=[tr tt];
B_270=[Br (BB*1e3)];
figure
hold on
grid on
plot(t_263,B_263,'r',t_270,B_270,'blue');
title(' Ice and water growth for Ta=263.15 K and Ta=270 K; Ts=Ta in both cases ');
legend('B@ 263.15K','B@ 270K','Position',[0.2, 0.78, .1, .1])
xlabel(' Time(s) '); ylabel(' B,h (mm) ');
  4 Comments
Adam Danz
Adam Danz on 27 Dec 2021
Upon a quick glance....
1. This line produces the red and blue lines. plot(t_263,B_263,'r',t_270,B_270,'blue');
2. Let's start with B_263 which is the red line. B_263 is defined by B_263 = [Br (B*1e3)];
3. If you plot Br and B*1e3 separately, you'll see that they form straight lines. So you are just combining two straight lines to form the red line. It's probably the same story for the blue line which I haven't looked at. figure(); plot(Br); figure(); plot(B*1e3).
My guess is that you're not using the same equations as the paper is using. Section VI of the paper points the equations used to create the ice and water thickness lines. I started mapping your code to the equations in the paper but realized it would take many more hours than I currently have time to invest. I suggest you work through the equations in the paper to understand which equations are used to plot fig 2 and after you fully understand that, start at the top of your code to determine where you're getting off track. This is where you become an expert in the matter 🤓.
mehmet salihi
mehmet salihi on 27 Dec 2021
you are right.
I thought about this method as I implement.
Even when i do not combine the two lines. the smooth line does not come with runga kutta method (Which has to come)

Sign in to comment.

Answers (0)

Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!