Can anybody help me to code boundary conditions in MATLAB for Keller Box Method?

Question

Prachi 2023-8-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2005302-can-anybody-help-me-to-code-boundary-conditions-in-matlab-for-keller-box-method

评论： Torsten 2024-2-27

f^'=1,f=S,θ^'=-r_1 [1+θ],ϕ^'=-r_2 [1+ϕ]               at            η=0
f^'=0,f^''=0,θ=0,ϕ=0                as                η→∞  

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mrutyunjaya Hiremath 2023-8-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2005302-can-anybody-help-me-to-code-boundary-conditions-in-matlab-for-keller-box-method#answer_1284002

在 MATLAB Online 中打开

% Define parameters

r_1 = 0.1;

r_2 = 0.2;

S = 2.0;

% Define the differential equations

% y(1) = f, y(2) = f', y(3) = θ, y(4) = ϕ

ode_system = @(eta, y) [y(2); 1; y(3); y(4)];

% Define the boundary conditions at η = 0

initial_conditions = [S, 1, 0, 0];

% Define the boundary conditions at η → ∞

eta_infinity = 100; % Choose a large value

final_conditions = [0, 0, 0, 0];

% Solve the differential equations

[eta, result] = ode45(ode_system, [0, eta_infinity], initial_conditions);

% Extract the solutions

f = result(:, 1);

f_prime = result(:, 2);

theta = result(:, 3);

phi = result(:, 4);

% Plot the solutions

subplot(2, 2, 1);

plot(eta, f);

xlabel('η');

ylabel('f');

title('f vs. η');

subplot(2, 2, 2);

plot(eta, f_prime);

xlabel('η');

ylabel("f'");

title("f' vs. η");

subplot(2, 2, 3);

plot(eta, theta);

xlabel('η');

ylabel('θ');

title('θ vs. η');

subplot(2, 2, 4);

plot(eta, phi);

xlabel('η');

ylabel('ϕ');

title('ϕ vs. η');

6 个评论
显示 4更早的评论隐藏 4更早的评论

Prachi 2023-8-7

I want to apply Keller box scheme for Jeffrey fluid.

Torsten 2023-8-7

It should be clear that we won't program this for you.

If you have a boundary value problem as above, you can use the MATLAB tools "bvp4c" or "bvp5c".

If your problem is an assignment, you will have to start programming it in MATLAB or make a google search whether you find a MATLAB code that fits your needs.

请先登录，再进行评论。

Answer 2

Santosh Devi 2024-2-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2005302-can-anybody-help-me-to-code-boundary-conditions-in-matlab-for-keller-box-method#answer_1417521

f^''' (η)+ff^'' (η)-(f^' (η))^2+Mf^' (η)-λf(η)=0

■θ^'' (η)+Pr⁡[f(η)θ^' (η)-b/(u_w^2 ) ηθ(η)+Ec(f^'' (η))^2+Q_0 θ(η)]=0

■ϕ^'' (η)+Sc[f(η) ϕ^' (η)-Kϕ(η)]=0

■f(0)=s,f^' (0)=1,θ(0)=1,ϕ(0)=1

■f^' (∞)→0,θ(∞)→0,ϕ(∞)→0

3 个评论
显示 1更早的评论隐藏 1更早的评论

Santosh Devi 2024-2-27

在 MATLAB Online 中打开

function plot_velocity_vs_eta()
    % Parameters
    lambda = 0.1;
    Pr = 0.72;
    b = 1;
    uw = 1;
    Ec = 0.0001;
    Q0 = 1;
    Sc = 1;
    K = 1;
    s = 1;
    M_values = [0, 0.5, 1];
    % Initialize plot
    figure;
    hold on;
    % Iterate over each Prandtl number
    for i = 1:length(M_values)
        M = M_values(i);
        % Define the system of equations
        equations = @(eta, y) [
            y(2);                                                       % f'(eta) = y(2)
            y(3);                                                       % f''(eta) = y(3)
            -y(1)*y(3) + y(2)^2 - M*y(2) + lambda*y(1);                % theta'(eta) = y(3)
            -Pr * (y(1)*y(4) - b/(uw^2)*eta*y(3) + Ec*y(3)^2 + Q0*y(4)); % theta''(eta)
            y(6);                                                       % phi'(eta) = y(6)
            -Sc * (y(1)*y(6) - K*y(6))                                  % phi''(eta)
        ];
        % Define the boundary conditions function
        bc = @(ya, yb) boundary_conditions(ya, yb, s);
        % Solve the equations
        eta_span = [0, 10];
        initial_conditions = [s; 0; 1; 1; 0; 0]; % Initial guesses for f, f', theta, phi, f'', phi'
        options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
        [eta, solution] = ode45(equations, eta_span, initial_conditions, options);
        % Calculate velocity (f'(eta))
        velocity = solution(:, 2);
        % Plot velocity against eta
        plot(eta, velocity, '-', 'LineWidth', 2, 'DisplayName', sprintf('M = %.2f', M));
    end
    % Add labels and legend
    xlabel('\eta');
    ylabel('Velocity (f''(\eta))');
    title('Velocity Profile for Different M Numbers');
    legend('Location', 'best');
    grid on;
    hold off;
end
function res = boundary_conditions(ya, yb, s)
    % Boundary conditions
    res = [ya(1) - s; ya(2); ya(3) - 1; ya(5) - 1; ya(4); ya(6); yb(1) - s; yb(2); yb(3) - 1; yb(5) - 1; yb(4); yb(6)];
end