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)
Afficher commentaires plus anciens
Kushagra Saurabh
le 4 Sep 2024
Modifié(e) : Torsten
le 13 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 commentaires
Réponse acceptée
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
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
Plus de réponses (1)
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)
[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
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 ?
Voir également
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!