Filter löschen
Filter löschen

I'm trying to solve jeffry hamel equation using RK4 but it's not giving output as expected.

5 Ansichten (letzte 30 Tage)
%%%% 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.

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 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
with s = 0 and h(0) = -1.7 g(1) = -1.7 with s = 0 and h(0) = -2 g(1) = -2 with s = 6 and h(0) = -4.3 g(1) = -0.78216 with s = 6 and h(0) = -4.1 g(1) = -0.67057
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 Kommentare
Torsten
Torsten am 13 Sep. 2024
Bearbeitet: Torsten am 13 Sep. 2024
I can only find the boundary condition U + K*n*dU/dy = 0 at y = 1. So you want to set K = 0 ?

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Torsten
Torsten am 13 Sep. 2024
Bearbeitet: Torsten am 13 Sep. 2024
s = 6;
Kn = 0.1;
U20 = -1;
U2 = fsolve(@(U2)fun_shoot(U2,s,Kn),U20)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
U2 = -4.0971
[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 Kommentare
Kushagra Saurabh
Kushagra Saurabh am 13 Sep. 2024
Thank you so much Torsten, it would be helpful if you coud just help me know where I was wrong so that I can learn.
Thanks for your kind help.
Torsten
Torsten am 13 Sep. 2024
Bearbeitet: Torsten am 13 Sep. 2024
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 ?

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Structures finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by