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

3 vues (au cours des 30 derniers jours)
%%%% 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.

Réponse acceptée

Alan Stevens
Alan Stevens le 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 commentaires
Torsten
Torsten le 13 Sep 2024
Modifié(e) : Torsten le 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 ?

Connectez-vous pour commenter.

Plus de réponses (1)

Torsten
Torsten le 13 Sep 2024
Modifié(e) : Torsten le 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 commentaires
Kushagra Saurabh
Kushagra Saurabh le 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 le 13 Sep 2024
Modifié(e) : Torsten le 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 ?

Connectez-vous pour commenter.

Catégories

En savoir plus sur Structures dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by