Thermal BL Added

RK4_FalknerSkan.m → F1_RK4_FalknerSkan.m

@@ -1,6 +1,6 @@
 %% RK4 numerical solver modified for solving Falkner-Skan equation.
 
-function [f_next, g_next, h_next] = RK4_FalknerSkan(dh, d_eta, f, g, h, i)
+function [f_next, g_next, h_next] = F1_RK4_FalknerSkan(dh, d_eta, f, g, h, i)
 
     dg1 = d_eta * h(i);
     df1 = d_eta * g(i);

F2_RK4_ThermalBL.m

@@ -0,0 +1,26 @@
+%% RK4 numerical solver modified for solving Falkner-Skan equation.
+
+function [f_next, g_next] = F2_RK4_ThermalBL(dh, d_eta, f, g, i, Pr, f_blasius)
+
+    if i > 117
+        f_blasius(i) = f_blasius(117) + d_eta * (i - 117);
+    end
+
+    dg1 = d_eta * dh(Pr, f_blasius(i), g(i));
+    df1 = d_eta * g(i);
+
+    dg2 = d_eta * dh(Pr, f_blasius(i), g(i) + dg1 * 0.5);
+    df2 = d_eta * (g(i) + dg1 * 0.5);
+
+    dg3 = d_eta * dh(Pr, f_blasius(i), g(i) + dg2 * 0.5);
+    df3 = d_eta * (g(i) + dg2 * 0.5);
+
+    dg4 = d_eta * dh(Pr, f_blasius(i), g(i) + dg3 * 0.5);
+    df4 = d_eta * (g(i) + dg3);
+
+    g(i + 1) = g(i) + (dg1 + 2*dg2 + 2*dg3 + dg4) / 6;
+    f(i + 1) = f(i) + (df1 + 2*df2 + 2*df3 + df4) / 6;
+
+    f_next = f;
+    g_next = g;
+end

M1_FalknerSkan.m

@@ -3,17 +3,21 @@
 clear; clc;
 
 beta_arr = [-0.19; 0; 0.29; 1];
-d_eta = 0.001;
-eta_max = 15;
+d_eta = 0.05;
+eta_max = 10;
+Kp = 0.001;        % step size for Newton.
+acc_target = 1e-5; % target accuracy.
 
 eta_master = zeros(eta_max/d_eta, 4);
 f_master = eta_master;
 g_master = f_master;
 h_master = f_master;
+hh_master = f_master;
 h_inits = zeros(1,4);
 last_etas = zeros(1,4);
 
-for beta_index = 1:length(beta_arr)    
+
+for beta_index = 1:length(beta_arr)
 
     beta = beta_arr(beta_index);
 
@@ -23,6 +27,7 @@
     f = zeros(N,1); % f(eta)
     g = f;          % g(eta) = f'
     h = f;          % h(eta) = f''
+    hh = f;         % hh(eta) = f'''
 
     dh = @(f, g, h) - h * f - beta * (1 - g^2); % dh/d_eta obtained from f''' + ff'' + b*[1-f'^2]
 
@@ -31,15 +36,13 @@
     g_inf = 0;         % f'(inf) = 1 boundary value checker.
 
     iteration = 0;     % iteration count over initial guess.
-    Kp = 0.1;          % step size for Newton.
-    acc_target = 1e-3; % target accuracy.
 
     while true
         iteration = iteration + 1;
 
         for i = 1:(N - 1)
-            [f, g, h] = RK4_FalknerSkan(dh, d_eta, f, g, h, i);
-
+            [f, g, h] = F1_RK4_FalknerSkan(dh, d_eta, f, g, h, i);
+            hh(i) = (h(i + 1) - h(i)) / d_eta;
             if h(i) < acc_target
                 g_inf = g(i);   % set g_inf to the converged value.
                 c_index = i;
@@ -60,6 +63,7 @@
             f = f(1:c_index);
             g = g(1:c_index);
             h = h(1:c_index);
+            hh = hh(1:c_index);
             break
         end
     end
@@ -68,8 +72,15 @@
     f_master(1:c_index, beta_index) = f;
     g_master(1:c_index, beta_index) = g;
     h_master(1:c_index, beta_index) = h;
+    hh_master(1:c_index, beta_index) = hh;
     h_inits(beta_index) = h(1);
     last_etas(beta_index) = c_index;
 end
 
-run("M2_Visualizer.m")
+run("M2_Visualizer.m")
+
+f_blasius = f_master(1:last_etas(2),2);
+g_blasius = f_master(1:last_etas(2),2);
+h_blasius = f_master(1:last_etas(2),2);
+
+% save thermalBlasius.mat f_blasius g_blasius h_blasius

M2_Visualizer.m

@@ -1,7 +1,12 @@
 %% This .m file contains the visualization part of the HW and is seperated 
 % with the intention of decluttering the main code.
 
-vel_figure = figure(1);
+% Point of inflection identification
+
+[poi_val, poi_pos] = min(abs(hh_master(1:round(last_etas(1)*0.5),1)));
+
+%% f' plot
+figure(1);
 plot(eta_master(1:last_etas(1),1), g_master(1:last_etas(1),1),"Color", "#071013", "LineWidth", 1.5)
 hold on
 xlabel("\eta", 'FontSize', 10)
@@ -12,13 +17,21 @@
 plot(eta_master(1:last_etas(3),3), g_master(1:last_etas(3),3),"Color", "#B4DC7F", "LineWidth", 1.5)
 plot(eta_master(1:last_etas(4),4), g_master(1:last_etas(4),4),"Color", "#dc602e", "LineWidth", 1.5)
 
+%% plot the point of inflection
+
+plot(eta_master(poi_pos,1), g_master(poi_pos, 1) , "o","MarkerSize",7,"Color","red")
+quiver(eta_master(poi_pos,1) + 1, g_master(poi_pos, 1) - 0.1, -1, 0.1, "filled", "Color","k")
+
+legend("\beta = -0.19","\beta = 0", "\beta = 0.29", "\beta = 1", "Location","southeast")
+text(eta_master(poi_pos,1) + 1, g_master(poi_pos, 1) - 0.1, " Point of inflection @ \eta = " + eta_master(poi_pos,1))
 
-legend("\beta = 1","\beta = 0.29", "\beta = 0", "\beta = -0.19", "Location","southeast")
 ax = gca;
 ax.FontSize = 10;
 
+%% f'' plot
 figure(2)
 plot(eta_master(1:last_etas(1),1), h_master(1:last_etas(1),1),"Color", "#071013", "LineWidth", 1.5)
+
 hold on
 xlabel("\eta", 'FontSize', 10)
 ylabel("f''(\eta)", 'FontSize', 10)
@@ -27,19 +40,25 @@
 plot(eta_master(1:last_etas(2),2), h_master(1:last_etas(2),2),"Color", "#23b5d3", "LineWidth", 1.5)
 plot(eta_master(1:last_etas(3),3), h_master(1:last_etas(3),3),"Color", "#B4DC7F", "LineWidth", 1.5)
 plot(eta_master(1:last_etas(4),4), h_master(1:last_etas(4),4),"Color", "#dc602e", "LineWidth", 1.5)
-legend("\beta = 1","\beta = 0.29", "\beta = 0", "\beta = -0.19", "Location","northeast")
 
+legend("\beta = -0.19","\beta = 0", "\beta = 0.29", "\beta = 1", "Location","northeast")
+
+
+%% f''' plot
 figure(3)
-plot(eta_master(1:last_etas(1) - 1,1), diff(h_master(1:last_etas(1),1)),"Color", "#071013", "LineWidth", 1.5)
+plot(eta_master(1:last_etas(1),1), hh_master(1:last_etas(1),1),"Color", "#071013", "LineWidth", 1.5)
 hold on
 xlabel("\eta", 'FontSize', 10)
 ylabel("f'''(\eta)", 'FontSize', 10)
 grid on
 
-plot(eta_master(1:last_etas(2)-1,2), diff(h_master(1:last_etas(2),2)),"Color", "#23b5d3", "LineWidth", 1.5)
-plot(eta_master(1:last_etas(3)-1,3), diff(h_master(1:last_etas(3),3)),"Color", "#B4DC7F", "LineWidth", 1.5)
-plot(eta_master(1:last_etas(4)-1,4), diff(h_master(1:last_etas(4),4)),"Color", "#dc602e", "LineWidth", 1.5)
-legend("\beta = 1","\beta = 0.29", "\beta = 0", "\beta = -0.19", "Location","northeast")
-
+plot( eta_master(1:last_etas(2),2), hh_master(1:last_etas(2),2),"Color", "#23b5d3", "LineWidth", 1.5)
+plot( eta_master(1:last_etas(3),3), hh_master(1:last_etas(3),3),"Color", "#B4DC7F", "LineWidth", 1.5)
+plot( eta_master(1:last_etas(4),4), hh_master(1:last_etas(4),4),"Color", "#dc602e", "LineWidth", 1.5)
 
+ax = gca;
+ax.FontSize = 10;
 
+yline(0,"LineWidth",1);
+xline(eta_master(poi_pos,1),"k--");
+legend("\beta = -0.19","\beta = 0", "\beta = 0.29", "\beta = 1", "Location","southeast")

M3_ThermalBL.m

@@ -0,0 +1,80 @@
+clear; clc;
+
+load thermalBlasius.mat
+
+prandtl_arr = [0.02 0.9 1 90];
+
+d_eta = 0.05;
+eta_max = 30;
+Kp = 0.001;        % step size for Newton.
+acc_target = 1e-4; % target accuracy.
+
+
+eta_master = zeros(ceil(eta_max/d_eta), 4);
+t_master = eta_master;
+dt_master = t_master;
+
+
+dt_inits = zeros(1,4);
+last_etas = zeros(1,4);
+dh = @(p, f, dt) -p * f * dt;
+
+for pr_index = 1:length(prandtl_arr)
+    pr = prandtl_arr(pr_index);
+
+    eta = 0:d_eta:eta_max;
+    N = length(eta);
+
+    t = zeros(N,1); % t(eta) = θ
+    dt = t;         % dt(eta) = θ'
+
+    t(1) = 1; % θ(0) = 1
+
+    dt(1) = -0.2;       % initial guess.
+    dt_last = dt(1);
+    t_inf = 0;         % free-stream boundary condition
+
+    iteration = 0;     % iteration count over initial guess.
+
+    while true
+        iteration = iteration + 1;
+
+        for i = 1:(N - 1)
+            [t, dt] = F2_RK4_ThermalBL(dh, d_eta, t, dt, i, pr, f_blasius);
+            if abs(dt(i)) < acc_target
+                t_inf = t(i);
+                c_index = i;
+                break
+            end
+            if i == (N - 1)
+                t_inf = t(N - 1);
+                c_index = N - 1;
+            end
+        end
+
+
+        if abs(t_inf) > acc_target
+            t = zeros(N,1); % reset matrices for next iteration
+            t(1) = 1;
+            dt = t;
+            dt(1) = dt_last - t_inf * Kp;
+            dt_last = dt(1);
+        else
+            disp("Target accuracy reached. dt(0) = " + string(dt(1)))
+            disp(iteration + " iterations completed.")
+            eta = eta(1:c_index);
+            t = t(1:c_index);
+            dt = dt(1:c_index);
+            break
+        end
+    end
+
+    eta_master(1:c_index,pr_index) = eta;
+    t_master(1:c_index, pr_index) = t;
+    dt_master(1:c_index, pr_index) = dt;
+    dt_inits(pr_index) = dt(1);
+    last_etas(pr_index) = c_index;
+
+end
+
+run("M4_ThermalVisualizer.m")

M4_ThermalVisualizer.m

@@ -0,0 +1,37 @@
+%% This .m file contains the visualization part of the HW and is seperated 
+% with the intention of decluttering the main code.
+
+
+%% theta plot
+figure(1);
+plot(eta_master(1:last_etas(1),1), t_master(1:last_etas(1),1),"Color", "#071013", "LineWidth", 1.5)
+hold on
+xlabel("\eta", 'FontSize', 10)
+ylabel("\theta", 'FontSize', 10)
+grid on
+
+plot(eta_master(1:last_etas(2),2), t_master(1:last_etas(2),2),"Color", "#23b5d3", "LineWidth", 1.5)
+plot(eta_master(1:last_etas(3),3), t_master(1:last_etas(3),3),"Color", "#B4DC7F", "LineWidth", 1.5)
+plot(eta_master(1:last_etas(4),4), t_master(1:last_etas(4),4),"Color", "#dc602e", "LineWidth", 1.5)
+
+legend("Pr = 0.02","Pr = 0.9", "Pr = 1", "Pr = 90", "Location","northeast")
+
+ax = gca;
+ax.FontSize = 10;
+
+%% f'' plot
+figure(2)
+plot(eta_master(1:last_etas(1),1), dt_master(1:last_etas(1),1),"Color", "#071013", "LineWidth", 1.5)
+
+hold on
+xlabel("\eta", 'FontSize', 10)
+ylabel("f''(\eta)", 'FontSize', 10)
+grid on
+
+plot(eta_master(1:last_etas(2),2), dt_master(1:last_etas(2),2),"Color", "#23b5d3", "LineWidth", 1.5)
+plot(eta_master(1:last_etas(3),3), dt_master(1:last_etas(3),3),"Color", "#B4DC7F", "LineWidth", 1.5)
+plot(eta_master(1:last_etas(4),4), dt_master(1:last_etas(4),4),"Color", "#dc602e", "LineWidth", 1.5)
+
+legend("Pr = 0.02","Pr = 0.9", "Pr = 1", "Pr = 90", "Location","southeast")
+
+