why I do not receive smooth transition and smooth lines

Question

mehmet salihi 2021-12-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1618165-why-i-do-not-receive-smooth-transition-and-smooth-lines

评论： mehmet salihi 2021-12-27

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Adam Danz 2021-12-27

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 2021-12-27

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)

请先登录，再进行评论。

请先登录，再回答此问题。

why I do not receive smooth transition and smooth lines

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

why I do not receive smooth transition and smooth lines

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论