|
[Sponsors] |
![]() |
![]() |
#1 |
Member
Kim jeamin
Join Date: Jul 2023
Posts: 36
Rep Power: 3 ![]() |
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, ![]() ![]() v(:,1) = 0; v(:,Ny) = 0; % 좌우 경계 v(1, ![]() ![]() % 창문처럼 처리할 절점 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'); |
|
![]() |
![]() |
![]() |
![]() |
#2 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#3 | |
Member
Kim jeamin
Join Date: Jul 2023
Posts: 36
Rep Power: 3 ![]() |
Quote:
I will try to do as you said. Thank you for your response. ![]() |
||
![]() |
![]() |
![]() |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
How to make CFD code architechture flexible enough for future modifications? | aerosayan | Main CFD Forum | 6 | September 3, 2021 04:35 |
How to write a good cfd code in fortran ? | Dhairya | Main CFD Forum | 13 | May 26, 2017 04:29 |
From 2D compressible code to 3D code | David Liu | Main CFD Forum | 22 | June 26, 2012 17:59 |
user friendly cfd code | waqar | Main CFD Forum | 19 | August 18, 2000 16:31 |
Open source CFD code development, possible? | Dr. Yazid Bindar | Main CFD Forum | 27 | July 18, 2000 00:18 |