why I do not receive smooth transition and smooth lines

hello dear friends.

could you please explain me why i do not receive smooth transition like below graph

the graph is taken from

https://arc.aiaa.org/doi/10.2514/2.1312

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

[Adam Danz]

I've edited the question to run the code and produce the outputs.

[mehmet salihi]

thank you for your friendly help.

[Adam Danz]

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]

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)