Errors trying to get the right plots

1 次查看（过去 30 天）

Cesar Cardenas 2023-3-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1937474-errors-trying-to-get-the-right-plots

评论： VBBV 2023-3-30

I'm trying to get the plots shown in the enclosed image. Not sure how to fix the code. Attached the .txt data file. Any help will be appreciated. Thanks

close all

clear all

clc

% Inputs

A = 10;

m = 0.06;

L = 1;

nu = 2e-5; % Kinematic viscosity of the fluid

rho = 10; % Fluid density

x_start = 0.2;

x_end = 1;

dx = 0.001;

y_start = 0;

VelN = load('posm_profiles.txt');

y_high = round(max(VelN(:,1)),2);

y_end = y_high;

dy = 0.0001;

N = (x_end - x_start) / dx + 1

N = 801

M = (y_end - y_start) / dy + 1

M = 501

Re_x = A * L / nu;

% Initialize arrays

U = zeros(M, N);

Ue = zeros(1, N);

Theta = zeros(1, N);

cf = zeros(1, N);

y = linspace(y_start, y_end, M);

x = linspace(x_start, x_end, N);

%Inital Data Interpolation

Uef = @(XX)A*(XX/L).^m;

u = interp1(VelN(:,1),VelN(:,2),y)*Uef(x_start);

% Boundary conditions

U(1, :) = 0;

U(M, :) = 1;

% Solve for velocity profile using an explicit finite method

for i = 1:N

Ue(i) = (x(i) / L)^m * A;

theta = zeros(1, M);

for j = 2:M-1

F = [U(j-1, i), U(j, i), U(j+1, i)];

Uy = (F(3) - F(1)) / (2 * dy);

Uyy = (F(3) - 2 * F(2) + F(1)) / dy^2;

theta(j) = theta(j-1) + U(j, i) * dy / Ue(i);

cf(i) = 2 * nu * Uy / Ue(i)^2;

U(j, i+1) = U(j, i) + (U(j, i) * (F(3) - F(1)) - U(j, i)^2 * theta(j) + nu * Uyy) * dx / Ue(i)^2;

end

Theta(i) = trapz(y, 1 - U(:, i));

end

% Plot results

figure

plot(x, Theta,'o-')

title('Momentum thickness as a function of x')

xlabel('x')

ylabel('\theta')

figure

semilogy(x, cf,'o-')

title('Skin friction coefficient as a function of x')

xlabel('x')

ylabel('C_f')

% figure

% plot(y,u)

% title('Velocity Profile')

% xlabel('u [m/s]')

% ylabel('y [mm]')

%axis([-.001 .05 0 y_end])

figure

for i = 1:100:N

plot(U(:, i), y,'-')

hold on

end

title('Velocity profiles')

xlabel('u [m/s]')

ylabel('y [mm]')

%axis([-.001 .05 0 y_end])

legend('x = 0.2 m', 'x = 0.3 m', 'x = 0.4 m', 'x = 0.5 m', 'x = 0.6 m', 'x = 0.7 m', 'x = 0.8 m', 'x = 0.9 m', 'x = 1.0 m')

figure

for i = 1:100:N

U_over_Ue = U(:, i) / Ue(i);

y_over_theta = y / Theta(i);

plot(U_over_Ue, y_over_theta,'o-')

hold on

end

title('Velocity profiles normalized by momentum thickness')

xlabel('u/U_e')

ylabel('y/\theta')

axis([-0.001 .125 0 1.01])

legend('x = 0.2 m', 'x = 0.3 m', 'x = 0.4 m', 'x = 0.5 m', 'x = 0.6 m', 'x = 0.7 m', 'x = 0.8 m', 'x = 0.9 m', 'x = 1.0 m')

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

Cesar Cardenas 2023-3-29

Ok thank you, however, I have not been able to get the plots as shown in the image attached, so not sure how to fix it.?

VBBV 2023-3-30

在 MATLAB Online 中打开

When vector U are all zeros, its less likely, that velocity profile can be evaluated correctly,

There needs to be a condition where the vector U needs to be updated with right values

theta(j) = theta(j-1) + U(j, i) * dy / Ue(i);
% ------------------->>> theta vector are all zeros, conisder u(j) or some
% input taken from VelN variable or combination of both u(j) & VelN
% Recheck whether this equation is correctly formulated
U(j, i+1) = U(j, i) + (U(j, i) * (F(3) - F(1)) - U(j, i)^2 * theta(j) + nu * Uyy) * dx / Ue(i)^2;