[mongya]

I am trying to code simple algorithm to solve 2 dimensional convection-diffusion problem(steady state). and what i wanted and expected was velocity vector and heat visualization assemble with driven cavity flow. but what the results are the pictures below. picture with vector is graph of velocity and the coloured one is graph of heat. the effects of convection doesn not seem to be reflected in the pictures at all. IMG_0845.jpgIMG_0846.jpg

The code :
% 초기 조건 설정
Nx = 7;
Ny = 9;
dx = 1.0;
dy = 1.0;
rho = 1.293; % kg/m^3
mu = 1.81e-5; % kg/(m·s)
k = 0.0262; % 열전도도 (W/m·K)
Cp = 1005; % 비열 (J/kg·K)
alpha = k / (rho * Cp); % 열 확산 계수 (m^2/s)
max_iter = 1000;
tolerance = 1e-7;
T_initial_C = 28; % 초기 온도 (섭씨)
T_initial = T_initial_C + 273.15; % 켈빈 온도로 변환
u_inlet = 5.0; % 경계 속도 (m/s)
v_inlet = 0.0; % y-방향 속도는 0으로 설정
alpha_u = 0.7; % 속도 이완 계수
alpha_p = 0.3; % 압력 이완 계수

% 그리드 생성
p = ones(Nx, Ny) * 101325; % 초기 압력 (1 atm)
u = zeros(Nx, Ny); % x-방향 속도
v = zeros(Nx, Ny); % y-방향 속도
T = ones(Nx, Ny) * T_initial; % 초기 온도 분포

% inlet 및 outlet 절점 설정
u(1, 5) = u_inlet; % (1,5) 절점에서 속도 설정 (inlet)
u(7, 5) = -u_inlet; % (7,5) 절점에서 속도 설정 (outlet)
v(1, 5) = v_inlet;
v(7, 5) = v_inlet;

% 나머지 경계 조건을 벽으로 설정 (속도 = 0)
u(:,1) = 0; u(:,Ny) = 0; % 좌우 경계
u(1, = 0; u(Nx, = 0; % 상하 경계
v(:,1) = 0; v(:,Ny) = 0; % 좌우 경계
v(1, = 0; v(Nx, = 0; % 상하 경계

% 창문처럼 처리할 절점
windows = [1,2; 1,3; 1,7; 1,8; 7,2; 7,3; 7,7; 7,8];

% 창문 절점에서의 속도 설정 (이전 값을 유지)
for k = 1:size(windows, 1)
i = windows(k, 1);
j = windows(k, 2);
u(i, j) = 0; % 창문에서 속도 초기화
v(i, j) = 0;
end

% 그리드 한 가운데 절점에 고정 온도 291K 설정
center_i = floor(Nx / 2) + 1;
center_j = floor(Ny / 2) + 1;
T(center_i, center_j) = 291; % 고정 온도 (K)

% SIMPLE 알고리즘 반복
for iter = 1:max_iter
% STEP 1: 모멘텀 방정식 해결
u_star = u;
v_star = v;
for i = 2:Nx-1
for j = 2:Ny-1
u_star(i, j) = u(i, j) + ...
(- (p(i+1, j) - p(i-1, j)) / (2 * rho * dx) + ...
mu * ((u(i+1, j) - 2*u(i, j) + u(i-1, j)) / dx^2 + ...
(u(i, j+1) - 2*u(i, j) + u(i, j-1)) / dy^2) / rho);
v_star(i, j) = v(i, j) + ...
(- (p(i, j+1) - p(i, j-1)) / (2 * rho * dy) + ...
mu * ((v(i+1, j) - 2*v(i, j) + v(i-1, j)) / dx^2 + ...
(v(i, j+1) - 2*v(i, j) + v(i, j-1)) / dy^2) / rho);
end
end

% STEP 2: 압력 보정 방정식 해결
p_prime = zeros(Nx, Ny);
for i = 2:Nx-1
for j = 2:Ny-1
p_prime(i, j) = (1 / (2 * (dx^2 + dy^2))) * ...
(rho * (dx^2 * dy^2) * ...
((u_star(i+1, j) - u_star(i-1, j)) / (2 * dx) + ...
(v_star(i, j+1) - v_star(i, j-1)) / (2 * dy)));
end
end

% STEP 3: 속도 및 압력 보정 (이완 계수 적용)
for i = 2:Nx-1
for j = 2:Ny-1
u(i, j) = u_star(i, j) - alpha_u * (p_prime(i+1, j) - p_prime(i-1, j)) / (2 * rho * dx);
v(i, j) = v_star(i, j) - alpha_u * (p_prime(i, j+1) - p_prime(i, j-1)) / (2 * rho * dy);
end
end

% 경계 조건 재설정 (속도 경계 조건 유지)
u(1, 5) = u_inlet; % (1,5) 절점에서 속도 설정 (inlet)
u(7, 5) = -u_inlet; % (7,5) 절점에서 속도 설정 (outlet)
v(1, 5) = v_inlet;
v(7, 5) = v_inlet;

% 압력 업데이트 (이완 계수 적용)
p = p + alpha_p * p_prime;

% 수렴 검사
if max(max(abs(p_prime))) < tolerance
fprintf('수렴 도달: %d 반복 후\n', iter);
break;
end
end

% STEP 4: 열 유속 계산 및 온도 분포 업데이트 (정상 상태)
for iter = 1:max_iter
T_new = T;
for i = 2:Nx-1
for j = 2:Ny-1
if i == center_i && j == center_j
continue; % 중앙 고정 온도 유지
end
T_new(i, j) = T(i, j) + ...
(alpha * ((T(i+1, j) - 2*T(i, j) + T(i-1, j)) / dx^2 + ...
(T(i, j+1) - 2*T(i, j) + T(i, j-1)) / dy^2) - ...
(u(i, j) * (T(i+1, j) - T(i-1, j)) / (2 * dx)) - ...
(v(i, j) * (T(i, j+1) - T(i, j-1)) / (2 * dy)));
end
end
% 수렴 검사
if max(max(abs(T_new - T))) < tolerance
fprintf('온도 수렴 도달: %d 반복 후\n', iter);
break;
end
T = T_new; % Update temperature for the next iteration
end

% 온도를 섭씨로 변환
T_C = T - 273.15;

% 결과 시각화
[X, Y] = meshgrid(1:Nx, 1:Ny);

% 온도 분포 시각화
figure;
contourf(X, Y, T_C', 20);
title('온도 분포 (섭씨)');
xlabel('X');
ylabel('Y');
colorbar;

% 속도 벡터 시각화
figure;
quiver(X, Y, u', v');
title('속도 벡터 분포');
xlabel('X');
ylabel('Y');

[FMDenaro]

I think you are too overly ambitious ... and no one will check all the lines of your code.

Start first showing you are able to solve a steady linear 2D convection-diffusion problem with mixed Dirichlet-Neumann BCs.



However, your way of programming need to be reformulated. If you have a domain of leghts Lx, Ly then you first decide the number of steps along x and y, for example Nx and Ny. Then
dx= Lx/Nx and dy=Ly/Ny.


That works both in a dimensional and non-dimensional solvers. Your decision to set dx=dy=1 in dimensional form makes no sense.



When you demonstrate you are able to solve simple linear problem, maybe you can afford a simple solver for NSE.

[mongya]

Originally Posted by FMDenaro

I think you are too overly ambitious ... and no one will check all the lines of your code.

Start first showing you are able to solve a steady linear 2D convection-diffusion problem with mixed Dirichlet-Neumann BCs.



However, your way of programming need to be reformulated. If you have a domain of leghts Lx, Ly then you first decide the number of steps along x and y, for example Nx and Ny. Then
dx= Lx/Nx and dy=Ly/Ny.


That works both in a dimensional and non-dimensional solvers. Your decision to set dx=dy=1 in dimensional form makes no sense.



When you demonstrate you are able to solve simple linear problem, maybe you can afford a simple solver for NSE.

I am Sorry Mr. FMDenaro. I should have been more careful when posting such posts.

I will try to do as you said. Thank you for your response.