Navier-Stokes Solver With Finite Difference Method - Unexpected Results

28 次查看（过去 30 天）
AAm
AAm 2023-3-24
评论： ASMA 2024-5-11
Ran in:
Dear All,
I have been working on a Matlab script to implement a Navier-Stokes solver for the lid-driven wall problem but I think there is something wrong in my math because I am not getting the expected results if I compare my plots with benchmarks I have collected from online literatures and also I am not able to infer turbolence in my results by increasing the Re number.
I have also added a little function that checks for the boundaries conditions.
The solver is using the finite difference method to update the velocity and pressure fields, it calculates the tentative velocity using the momentum equations, corrects the pressure using the Poisson equation, finally it updates the velocity by using the corrected pressure. The boundary conditions are also applied.
This is the code if any of you fancy to help!
Alex
% Navier-Stokes Solver for Incompressible, Isothermal, and Newtonian Flow
clc
% Parameters
Lx = 1; % Length in x-direction
Ly = 1; % Length in y-direction
Nx = 32; % Number of grid points in x-direction
Ny = 32; % Number of grid points in y-direction
Re = 100; % Reynolds number
U_wall = 1; % Lid-driven wall velocity
tolerance = 1e-4; % Convergence tolerance
omega = 10; % frequency of oscillations
Amp = 0.1; % amplitude of oscillations
% Discretization
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Initialize velocity and pressure fields
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
p = zeros(Ny, Nx);
% Time stepping parameters
dt = 0.001;
max_iter = 50000;
% Main loop
for iter = 1:max_iter
% Calculate lid-driven wall velocity with sinusoidal variations
U_wall = 1 + Amp*sin(omega*dt*iter);
[u, v, p] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance);
if mod(iter, 100) == 0
fprintf('Iteration: %d\n', iter);
pause(0.01);
end
end
Iteration: 100 Iteration: 200 Iteration: 300 Iteration: 400 Iteration: 500 Iteration: 600 Iteration: 700 Iteration: 800 Iteration: 900 Iteration: 1000 Iteration: 1100 Iteration: 1200 Iteration: 1300 Iteration: 1400 Iteration: 1500 Iteration: 1600 Iteration: 1700 Iteration: 1800 Iteration: 1900 Iteration: 2000 Iteration: 2100 Iteration: 2200 Iteration: 2300 Iteration: 2400 Iteration: 2500 Iteration: 2600 Iteration: 2700 Iteration: 2800 Iteration: 2900 Iteration: 3000 Iteration: 3100 Iteration: 3200 Iteration: 3300 Iteration: 3400 Iteration: 3500 Iteration: 3600 Iteration: 3700 Iteration: 3800 Iteration: 3900 Iteration: 4000 Iteration: 4100 Iteration: 4200 Iteration: 4300 Iteration: 4400 Iteration: 4500 Iteration: 4600 Iteration: 4700 Iteration: 4800 Iteration: 4900 Iteration: 5000 Iteration: 5100 Iteration: 5200 Iteration: 5300 Iteration: 5400 Iteration: 5500 Iteration: 5600 Iteration: 5700 Iteration: 5800 Iteration: 5900 Iteration: 6000 Iteration: 6100 Iteration: 6200 Iteration: 6300 Iteration: 6400 Iteration: 6500 Iteration: 6600 Iteration: 6700 Iteration: 6800 Iteration: 6900 Iteration: 7000 Iteration: 7100 Iteration: 7200 Iteration: 7300 Iteration: 7400 Iteration: 7500 Iteration: 7600 Iteration: 7700 Iteration: 7800 Iteration: 7900 Iteration: 8000 Iteration: 8100 Iteration: 8200 Iteration: 8300 Iteration: 8400 Iteration: 8500 Iteration: 8600 Iteration: 8700 Iteration: 8800 Iteration: 8900 Iteration: 9000 Iteration: 9100 Iteration: 9200 Iteration: 9300 Iteration: 9400 Iteration: 9500 Iteration: 9600 Iteration: 9700 Iteration: 9800 Iteration: 9900 Iteration: 10000 Iteration: 10100 Iteration: 10200 Iteration: 10300 Iteration: 10400 Iteration: 10500 Iteration: 10600 Iteration: 10700 Iteration: 10800 Iteration: 10900 Iteration: 11000 Iteration: 11100 Iteration: 11200 Iteration: 11300 Iteration: 11400 Iteration: 11500 Iteration: 11600 Iteration: 11700 Iteration: 11800 Iteration: 11900 Iteration: 12000 Iteration: 12100 Iteration: 12200 Iteration: 12300 Iteration: 12400 Iteration: 12500 Iteration: 12600 Iteration: 12700 Iteration: 12800 Iteration: 12900 Iteration: 13000 Iteration: 13100 Iteration: 13200 Iteration: 13300 Iteration: 13400 Iteration: 13500 Iteration: 13600 Iteration: 13700 Iteration: 13800 Iteration: 13900 Iteration: 14000 Iteration: 14100 Iteration: 14200 Iteration: 14300 Iteration: 14400 Iteration: 14500 Iteration: 14600 Iteration: 14700 Iteration: 14800 Iteration: 14900 Iteration: 15000 Iteration: 15100 Iteration: 15200 Iteration: 15300 Iteration: 15400 Iteration: 15500 Iteration: 15600 Iteration: 15700 Iteration: 15800 Iteration: 15900 Iteration: 16000 Iteration: 16100 Iteration: 16200 Iteration: 16300 Iteration: 16400 Iteration: 16500 Iteration: 16600 Iteration: 16700 Iteration: 16800 Iteration: 16900 Iteration: 17000 Iteration: 17100 Iteration: 17200 Iteration: 17300 Iteration: 17400 Iteration: 17500 Iteration: 17600 Iteration: 17700 Iteration: 17800 Iteration: 17900 Iteration: 18000 Iteration: 18100 Iteration: 18200 Iteration: 18300 Iteration: 18400 Iteration: 18500 Iteration: 18600 Iteration: 18700 Iteration: 18800 Iteration: 18900 Iteration: 19000 Iteration: 19100 Iteration: 19200 Iteration: 19300 Iteration: 19400 Iteration: 19500 Iteration: 19600 Iteration: 19700 Iteration: 19800 Iteration: 19900 Iteration: 20000 Iteration: 20100 Iteration: 20200 Iteration: 20300 Iteration: 20400 Iteration: 20500 Iteration: 20600 Iteration: 20700 Iteration: 20800 Iteration: 20900 Iteration: 21000 Iteration: 21100 Iteration: 21200 Iteration: 21300 Iteration: 21400 Iteration: 21500 Iteration: 21600 Iteration: 21700 Iteration: 21800 Iteration: 21900 Iteration: 22000 Iteration: 22100 Iteration: 22200 Iteration: 22300 Iteration: 22400 Iteration: 22500 Iteration: 22600 Iteration: 22700 Iteration: 22800 Iteration: 22900 Iteration: 23000 Iteration: 23100 Iteration: 23200 Iteration: 23300 Iteration: 23400 Iteration: 23500 Iteration: 23600 Iteration: 23700 Iteration: 23800 Iteration: 23900 Iteration: 24000 Iteration: 24100 Iteration: 24200 Iteration: 24300 Iteration: 24400 Iteration: 24500 Iteration: 24600 Iteration: 24700 Iteration: 24800 Iteration: 24900 Iteration: 25000 Iteration: 25100 Iteration: 25200 Iteration: 25300 Iteration: 25400 Iteration: 25500 Iteration: 25600 Iteration: 25700 Iteration: 25800 Iteration: 25900 Iteration: 26000 Iteration: 26100 Iteration: 26200 Iteration: 26300 Iteration: 26400 Iteration: 26500 Iteration: 26600 Iteration: 26700 Iteration: 26800 Iteration: 26900 Iteration: 27000 Iteration: 27100 Iteration: 27200 Iteration: 27300 Iteration: 27400 Iteration: 27500 Iteration: 27600 Iteration: 27700 Iteration: 27800 Iteration: 27900 Iteration: 28000 Iteration: 28100 Iteration: 28200 Iteration: 28300 Iteration: 28400 Iteration: 28500 Iteration: 28600 Iteration: 28700 Iteration: 28800 Iteration: 28900 Iteration: 29000 Iteration: 29100 Iteration: 29200 Iteration: 29300 Iteration: 29400 Iteration: 29500 Iteration: 29600 Iteration: 29700 Iteration: 29800 Iteration: 29900 Iteration: 30000 Iteration: 30100 Iteration: 30200 Iteration: 30300 Iteration: 30400 Iteration: 30500 Iteration: 30600 Iteration: 30700 Iteration: 30800 Iteration: 30900 Iteration: 31000 Iteration: 31100 Iteration: 31200 Iteration: 31300 Iteration: 31400 Iteration: 31500 Iteration: 31600 Iteration: 31700 Iteration: 31800 Iteration: 31900 Iteration: 32000 Iteration: 32100 Iteration: 32200 Iteration: 32300 Iteration: 32400 Iteration: 32500 Iteration: 32600 Iteration: 32700 Iteration: 32800 Iteration: 32900 Iteration: 33000 Iteration: 33100 Iteration: 33200 Iteration: 33300 Iteration: 33400 Iteration: 33500 Iteration: 33600 Iteration: 33700 Iteration: 33800 Iteration: 33900 Iteration: 34000 Iteration: 34100 Iteration: 34200 Iteration: 34300 Iteration: 34400 Iteration: 34500 Iteration: 34600 Iteration: 34700 Iteration: 34800 Iteration: 34900 Iteration: 35000 Iteration: 35100 Iteration: 35200 Iteration: 35300 Iteration: 35400 Iteration: 35500 Iteration: 35600 Iteration: 35700 Iteration: 35800 Iteration: 35900 Iteration: 36000 Iteration: 36100 Iteration: 36200 Iteration: 36300 Iteration: 36400 Iteration: 36500 Iteration: 36600 Iteration: 36700 Iteration: 36800 Iteration: 36900 Iteration: 37000 Iteration: 37100 Iteration: 37200 Iteration: 37300 Iteration: 37400 Iteration: 37500 Iteration: 37600 Iteration: 37700 Iteration: 37800 Iteration: 37900 Iteration: 38000 Iteration: 38100 Iteration: 38200 Iteration: 38300 Iteration: 38400 Iteration: 38500 Iteration: 38600 Iteration: 38700 Iteration: 38800 Iteration: 38900 Iteration: 39000 Iteration: 39100 Iteration: 39200 Iteration: 39300 Iteration: 39400 Iteration: 39500 Iteration: 39600 Iteration: 39700 Iteration: 39800 Iteration: 39900 Iteration: 40000 Iteration: 40100 Iteration: 40200 Iteration: 40300 Iteration: 40400 Iteration: 40500 Iteration: 40600 Iteration: 40700 Iteration: 40800 Iteration: 40900 Iteration: 41000 Iteration: 41100 Iteration: 41200 Iteration: 41300 Iteration: 41400 Iteration: 41500 Iteration: 41600 Iteration: 41700 Iteration: 41800 Iteration: 41900 Iteration: 42000 Iteration: 42100 Iteration: 42200 Iteration: 42300 Iteration: 42400 Iteration: 42500 Iteration: 42600 Iteration: 42700 Iteration: 42800 Iteration: 42900 Iteration: 43000 Iteration: 43100 Iteration: 43200 Iteration: 43300 Iteration: 43400 Iteration: 43500 Iteration: 43600 Iteration: 43700 Iteration: 43800 Iteration: 43900 Iteration: 44000 Iteration: 44100 Iteration: 44200 Iteration: 44300 Iteration: 44400 Iteration: 44500 Iteration: 44600 Iteration: 44700 Iteration: 44800 Iteration: 44900 Iteration: 45000 Iteration: 45100 Iteration: 45200 Iteration: 45300 Iteration: 45400 Iteration: 45500 Iteration: 45600 Iteration: 45700 Iteration: 45800 Iteration: 45900 Iteration: 46000 Iteration: 46100 Iteration: 46200 Iteration: 46300 Iteration: 46400 Iteration: 46500 Iteration: 46600 Iteration: 46700 Iteration: 46800 Iteration: 46900 Iteration: 47000 Iteration: 47100 Iteration: 47200 Iteration: 47300 Iteration: 47400 Iteration: 47500 Iteration: 47600 Iteration: 47700 Iteration: 47800 Iteration: 47900 Iteration: 48000 Iteration: 48100 Iteration: 48200 Iteration: 48300 Iteration: 48400 Iteration: 48500 Iteration: 48600 Iteration: 48700 Iteration: 48800 Iteration: 48900 Iteration: 49000 Iteration: 49100 Iteration: 49200 Iteration: 49300 Iteration: 49400 Iteration: 49500 Iteration: 49600 Iteration: 49700 Iteration: 49800 Iteration: 49900 Iteration: 50000
%Plot Pressure Field
contourf(X, Y, p);
title('Pressure');
colorbar;
pause(0.01);
% Plot velocity field
figure;
quiver(X, Y, u, v);
title('Velocity Field');
xlabel('x');
ylabel('y');
figure;
contourf(X,Y,u',23,'LineColor','none');
title('U-Velocity');
colorbar;
%Plot stream function
N = 32; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
%[X,Y] = meshgrid(x,y);
figure; h=streamline(X,Y,u',v',xstart,ystart,[0.1, 200]);
title('Stream Function'); xlabel('x-location'); ylabel('y-location')
axis('equal',[0 Lx 0 Ly]); set(h,'color','k')
boundary_condition_tests(u, v, U_wall);
Boundary Condition Test Bottom wall u-velocity: Fail (Actual value: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.9064456389355]) Top wall u-velocity: Fail (Actual value: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.9064456389355]) Left and right wall u-velocity: Fail (Actual value: Left - [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0], Right - [1.9064456389355;1.26887379000312;1.07979496727235;1.01050634137285;0.979186948571896;0.962431256071441;0.95235177731503;0.945813692801267;0.94137177247205;0.93827141994475;0.936078123607436;0.934523755206433;0.933435685907068;0.932701026704101;0.932247075247037;0.932030307327377;0.93203030732738;0.932247075247041;0.932701026704096;0.933435685907069;0.934523755206432;0.936078123607436;0.938271419944747;0.941371772472047;0.94581369280126;0.952351777315031;0.962431256071432;0.979186948571894;1.01050634137285;1.07979496727235;1.26887379000311;1.9064456389355]) Left and right wall v-velocity: Pass Bottom and top wall v-velocity: Pass Total failed boundary conditions: 3
function [u_new, v_new, p_new] = navier_stokes_solver(u, v, p, dx, dy, dt, Re, U_wall, tolerance)
[Ny, Nx] = size(u);
% Initialize temporary fields
u_new = u;
v_new = v;
p_new = p;
% Calculate the tentative velocity (u*) using momentum equations
for i = 2:Ny - 1
for j = 2:Nx - 1
% Viscous acceleration term (in x-direction)
viscous_x_term = (1 / Re) * (u(i, j + 1) - 2 * u(i, j) + u(i, j - 1)) / (dx^2);
% Viscous acceleration term (in y-direction)
viscous_y_term = (1 / Re) * (u(i + 1, j) - 2 * u(i, j) + u(i - 1, j)) / (dy^2);
u_new(i, j) = u(i, j) + dt * (viscous_x_term + viscous_y_term);
% Viscous term in x-direction
viscous_v_x_term = (1 / Re) * (v(i, j + 1) - 2 * v(i, j) + v(i, j - 1)) / (dx^2);
viscous_v_y_term = (1 / Re) * (v(i + 1, j) - 2 * v(i, j) + v(i - 1, j)) / (dy^2);
v_new(i, j) = v(i, j) + dt * (viscous_v_x_term + viscous_v_y_term);
end
end
% Pressure correction (Poisson equation)
p_initial = p_new; % Store the initial pressure field
maxe = 1;
while maxe > tolerance
p_old = p_new;
for i = 2:Ny - 1
for j = 2:Nx - 1
p_new(i, j) = (1 / 4) * (p_old(i, j + 1) + p_old(i, j - 1) ...
+ p_old(i + 1, j) + p_old(i - 1, j) ...
- dx * (u_new(i, j + 1) - u_new(i, j - 1)) / 2 ...
- dy * (v_new(i + 1, j) - v_new(i - 1, j)) / 2);
end
end
maxe = max(max(abs(p_new - p_old)));
end
%Update velocity with pressure correction using the initial pressure field
for i = 2:Ny - 1
for j = 2:Nx - 1
u_new(i, j) = u_new(i, j) - dt * (p_initial(i, j + 1) - p_initial(i, j - 1)) / (2 * dx);
v_new(i, j) = v_new(i, j) - dt * (p_initial(i + 1, j) - p_initial(i - 1, j)) / (2 * dy);
end
end
% Boundary conditions
u_new(1, :) = 0;
u_new(end, :) = 0;
u_new(:, 1) = 0;
u_new(:, end) = 2 * U_wall - u_new(:, end - 1); % Moving lid
v_new(1, :) = -v_new(2, :);
v_new(end, :) = -v_new(end - 1, :);
v_new(:, 1) = 0;
v_new(:, end) = 0;
p_new(1, :) = p_new(2, :);
p_new(end, :) = p_new(end - 1, :);
p_new(:, 1) = p_new(:, 2);
p_new(:, end) = p_new(:, end - 1);
end
function boundary_condition_tests(u, v, U_wall)
fprintf('Boundary Condition Test\n');
% Bottom wall u-velocity
bottom_wall_u = u(1, :);
bottom_wall_u_test = all(abs(bottom_wall_u) < 1e-6);
if bottom_wall_u_test
fprintf('Bottom wall u-velocity: Pass\n');
else
fprintf('Bottom wall u-velocity: Fail (Actual value: %s)\n', mat2str(bottom_wall_u));
end
% Top wall u-velocity
top_wall_u = u(end, :);
top_wall_u_test = all(abs(top_wall_u - 2 * U_wall) < 1e-6);
if top_wall_u_test
fprintf('Top wall u-velocity: Pass\n');
else
fprintf('Top wall u-velocity: Fail (Actual value: %s)\n', mat2str(top_wall_u));
end
% Left and right wall u-velocity
left_wall_u = u(:, 1);
right_wall_u = u(:, end);
left_right_wall_u_test = all(abs(left_wall_u) < 1e-6) && all(abs(right_wall_u) < 1e-6);
if left_right_wall_u_test
fprintf('Left and right wall u-velocity: Pass\n');
else
fprintf('Left and right wall u-velocity: Fail (Actual value: Left - %s, Right - %s)\n', mat2str(left_wall_u), mat2str(right_wall_u));
end
% Left and right wall v-velocity
left_wall_v = v(:, 1);
right_wall_v = v(:, end);
left_right_wall_v_test = all(abs(left_wall_v) < 1e-6) && all(abs(right_wall_v) < 1e-6);
if left_right_wall_v_test
fprintf('Left and right wall v-velocity: Pass\n');
else
fprintf('Left and right wall v-velocity: Fail (Actual value: Left - %s, Right - %s)\n', mat2str(left_wall_v), mat2str(right_wall_v));
end
% Bottom and top wall v-velocity
bottom_wall_v = v(1, :);
top_wall_v = v(end, :);
bottom_top_wall_v_test = all(abs(bottom_wall_v + v(2, :)) < 1e-6) && all(abs(top_wall_v + v(end - 1, :)) < 1e-6);
if bottom_top_wall_v_test
fprintf('Bottom and top wall v-velocity: Pass\n');
else
fprintf('Bottom and top wall v-velocity: Fail (Actual value: Bottom - %s, Top - %s)\n', mat2str(bottom_wall_v), mat2str(top_wall_v));
end
% Count the total failed boundary conditions
failed_count = sum([~bottom_wall_u_test, ~top_wall_u_test, ~left_right_wall_u_test, ~left_right_wall_v_test, ~bottom_top_wall_v_test]);
fprintf('Total failed boundary conditions: %d\n', failed_count);
end

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

回答（1 个）

Alessio
Alessio 2024-2-1
In computing the u* velocity fields, you only have the diffusive/viscous term but the convective term is missing.

类别

在 Help Center 和 File Exchange 中查找有关 Oceanography and Hydrology 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by