I'm trying to solve jeffry hamel equation using RK4 but it's not giving output as expected.
44 views (last 30 days)
Show older comments
Kushagra Saurabh
on 4 Sep 2024
%%%% Jeffry-Hamel equation stability analysis %%%%
%%% ode:(U''' +2sUU' = 0) such that U' = g, g' = h, h' = -2*s*U*g %%%
%%% where s = S/a and a = slope of the wall %%%
%%% BC: U =1,g =0 at y =0 && U = 0, g= 10 (Kn = 0.1) at y = 1 %%%
%%% Find h(0) such that g(1) = 10 %%%
%% Range of Variable %%
y_cl = 0;
y_chw = 1;
Re = 0.68e5;
S = 6;
%% Boundary conditions %%
U0 = 1;
U1 = 0;
g0 = 0;
g1 = 10;
h0 = [0.1 0.3];
%% Improved step size and error tolerance %%
d_y = 0.001; % Reduce step size for better accuracy
N = (y_chw -y_cl)/d_y;
err = 1;
max_iter = 6e5; % Limit the maximum number of iterations
iter_count = 0;
%% initializing solution %%
y = y_cl : d_y : y_chw; % y coordinate for plot function
while err > 1e-15 && iter_count < max_iter
iter_count = iter_count + 1;
for j = 1:2
F = zeros(N+1, 3); % Initialize solution matrix
F(1, :) = [U0 g0 h0(j)];
% Runge-Kutta 4th order method
for i = 1:N
k1 = d_y * [F(i, 2) F(i, 3) -(F(i, 1) * F(i, 2)) * 2 * S];
k2 = d_y * [F(i, 2) + k1(2)/2 F(i, 3) + k1(3)/2 -( (F(i, 1) + k1(1)/2) * (F(i, 2) + k1(2)/2) ) * 2 * S];
k3 = d_y * [F(i, 2) + k2(2)/2 F(i, 3) + k2(3)/2 -( (F(i, 1) + k2(1)/2) * (F(i, 2) + k2(2)/2) ) * 2 * S];
k4 = d_y * [F(i, 2) + k3(2) F(i, 3) + k3(3) -( (F(i, 1) + k3(1)) * (F(i, 2) + k3(2)) ) * 2 * S];
F(i+1, :) = F(i, :) + (1/6) * (k1 + 2*k2 + 2*k3 + k4);
end
g_1(j) = F(end, 2); % Store g(1) for each trial
end
% Compute error and update h0
[err, index] = max(abs(g1 - g_1));
P = diff(g_1);
if P ~= 0
h0_new = h0(1) + (diff(h0) / diff(g_1)) * (g1 - g_1(1));
h0(index) = h0_new;
end
end
% Final check on iteration count
if iter_count >= max_iter
warning('Max iterations reached. Convergence may not be achieved.');
end
%% plotting
f1=figure;
f1.Units = 'normalized';
f1.Position = [0.1 0.1 0.8 0.6];
plot(F,y','LineWidth',1);
xlim([0 1]);
ylim([0 1]);
axis square
yticks(0:0.2:1);
yticklabels({'0','0.2','0.4','0.6','0.8','1'});
xticks(0:0.2:1);
xticklabels({'0','0.2','0.4','0.6','0.8','1'});
grid on
title('Jeffry Hamel solution');
xlabel('y');
legend('U','U"','U"''');
I'm also inserting a plot the I want to see:
please help me in getting this resolved.
Thanks in advance.
0 Comments
Accepted Answer
Alan Stevens
on 6 Sep 2024
The following code produces something very close to the plot that you want to see. However, the values of g(1) are nowhere near 10. Are you sure you have expressed the Jeffery-Hamel equations correctly?
%%%% Jeffry-Hamel equation stability analysis %%%%
%%% ode:(U''' +2sUU' = 0) such that U' = g, g' = h, h' = -2*s*U*g %%%
%%% where s = S/a and a = slope of the wall %%%
%%% BC: U =1,g =0 at y =0 && U = 0, g= 10 (Kn = 0.1) at y = 1 %%%
%%% Find h(0) such that g(1) = 10 %%%
yspan = 0:0.1:1;
U0 = 1; g0 = 0;
h0 = [-1.7, -2.0, -4.3, -4.1];
s = [0, 0, 6, 6];
sl = ['-s'; '-s'; '-+'; '-+'];
figure
hold on
for i = 1:numel(h0)
Ugh0 = [U0, g0, h0(i)];
[y, Ugh] = ode45(@(y,Ugh)fn(y,Ugh,s(i)), yspan, Ugh0);
U = Ugh(:,1); g = Ugh(:,2); h = Ugh(:,3);
plot(U,y,sl(i,:))
disp(['with s = ', num2str(s(i)),' and h(0) = ',num2str(h0(i)),...
' g(1) = ', num2str(g(end))])
end
grid
xlabel('U'), ylabel('y')
axis([0,1,0,1])
legend(['s = 0', ' h(0) = ', num2str(h0(1))], ...
['s = 0 ',' h(0) = ', num2str(h0(2))], ...
['s = 6', ' h(0) = ', num2str(h0(3))], ...
['s = 6', ' h(0) = ', num2str(h0(4))])
function dUghdy = fn(~,Ugh,s)
U = Ugh(1); g = Ugh(2); h = Ugh(3);
dUghdy = [g; h; -2*s*U*g];
end
11 Comments
More Answers (1)
Torsten
on 13 Sep 2024 at 10:18
Edited: Torsten
on 13 Sep 2024 at 10:29
s = 6;
Kn = 0.1;
U20 = -1;
U2 = fsolve(@(U2)fun_shoot(U2,s,Kn),U20)
[Y,Z] = fun_odes(U2,s);
plot(Z(:,1),Y)
xlabel('U')
ylabel('y')
grid on
function res = fun_shoot(U2,s,Kn)
[Y,Z] = fun_odes(U2,s);
res = Z(end,1) + Kn* Z(end,2);
end
function [Y,Z] = fun_odes(U2,s)
f = @(y,z)[z(2);z(3);-2*s*z(1)*z(2)];
z0 = [1;0;U2];
[Y,Z] = ode45(f,[0 1],z0);
end
2 Comments
Torsten
on 13 Sep 2024 at 10:48
Edited: Torsten
on 13 Sep 2024 at 12:24
I don't know - I didn't check your code. I just used MATLAB functions ("fsolve" for shooting and "ode45" for integration) and it worked fine.
So what you should learn is to use as many MATLAB functions and to implement as little self-made numerical parts as possible to solve a problem.
E.g. the update
F(3) = F(3) - error; % Update h(0) based on error Converges quicker here if error not multiplied by a small factor
doesn't make sense at all in your code. How should the values you get for U at y=1 influence the 2nd derivative of U at y=0 in this way ?
See Also
Categories
Find more on Structures 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!