solution for steady flow at Re = 10. in coursera i took mathematics for engineers.for that i need to solve an assignment. For this code i cant get correct stream & vorticity
    조회 수: 21 (최근 30일)
  
       이전 댓글 표시
    
flow_around_cylinder_steady
function [psi, omega] = flow_around_cylinder_steady
    Re=10; 
    %%%%% define the grid %%%%%
    n=101; m=101; % number of grid points
    N=n-1; M=m-1; % number of grid intervals
    h=pi/M; % grid spacing based on theta variable
    xi=(0:N)*h; theta=(0:M)*h; % xi and theta variables on the grid
    %%%%% Initialize the flow fields %%%%%
    psi=zeros(n,m);
    omega=zeros(n,m);
    psi(n,:)=0 % Free stream boundary condition (psi = 0 at the top edge)
    %%%%% Set relax params, tol, extra variables %%%%%
    r_psi=1.8; % Set the relaxation parameter here, psi equation
    r_omega=0.9; % Set the relaxation parameter here, omega equation
    delta=1.e-08; % error tolerance
    error=2*delta; % initialize error variable
    %%%%% Add any additional variable definitions here %%%%%
    ...
    ...
    %%%%% Main SOR Loop %%%%%
    while (error > delta)
        psi_old = psi; omega_old = omega;
        for i=2:n-1
            for j=2:m-1
                psi(i,j)=... % (1 - r_psi) * psi_old(i, j) + ...
                    r_psi * 0.5 * (psi_old(i+1, j) + psi_old(i-1, j) + ...
                    psi_old(i, j+1) + psi_old(i, j-1) - h^2 * omega(i, j));
            end
        end
        error_psi=max(abs(psi(:)-psi_old(:)));
        omega(1,:)=0 % Boundary condition at theta = 0
        for i=2:n-1
            for j=2:m-1
                omega(i,j)=... % (1 - r_omega) * omega_old(i, j) + ...
                    r_omega * ( (omega_old(i+1, j) + omega_old(i-1, j) + ...
                    omega_old(i, j+1) + omega_old(i, j-1) - ...
                    (4 - (2/h^2)) * omega_old(i, j)) / (4 - (2/h^2)) );
            end
        end
        error_omega=max(abs(omega(:)-omega_old(:)));
        error=max(error_psi, error_omega);
    end
    plot_Re10(psi);
end
Please refer to page 59 (64/69) of the PDF document for the 'plot_Re10' custom plotting function.
댓글 수: 12
  Torsten
      
      
 2024년 8월 8일
				And where do you take care of the inscribed cylinder in your code ? I only see a rectangular region in which you compute "psi" and "omega".
답변 (3개)
  LeoAiE
      
 2024년 8월 7일
        
      편집: Sam Chak
      
      
 2024년 8월 8일
  
      Hi there! 
This is not exact solution you looking for but more to help you think about the issue you are facing! 
[psi, omega] = flow_around_cylinder_steady
function [psi, omega] = flow_around_cylinder_steady
    Re = 10;
    % Define the grid
    n = 101; 
    m = 101; 
    N = n - 1; 
    M = m - 1; 
    h = 2 * pi / M; 
    xi = linspace(0, 1, N+1); 
    theta = linspace(0, 2*pi, M+1); 
    % Initialize the flow fields
    psi = zeros(n, m);
    omega = zeros(n, m);
    % Boundary conditions for psi (streamfunction)
    psi(n, :) = 0; % Top boundary (free stream)
    psi(1, :) = 0; % Bottom boundary (cylinder surface)
    psi(:, 1) = psi(:, 2); % Left boundary (symmetry)
    psi(:, end) = psi(:, end-1); % Right boundary (symmetry)
    % Boundary conditions for omega (vorticity)
    omega(1, :) = -2 * psi(2, :) / h^2; % No-slip condition at cylinder surface
    omega(:, 1) = omega(:, 2); % Left boundary (symmetry)
    omega(:, end) = omega(:, end-1); % Right boundary (symmetry)
    % Set relax params, tol, extra variables
    r_psi = 1.8; 
    r_omega = 0.9; 
    delta = 1.e-08; 
    error = 2 * delta; 
    % Main SOR Loop
    while (error > delta)
        psi_old = psi; 
        omega_old = omega;
        % Update psi
        for i = 2:n-1
            for j = 2:m-1
                psi(i, j) = (1 - r_psi) * psi_old(i, j) + ...
                    r_psi * 0.25 * (psi_old(i+1, j) + psi_old(i-1, j) + ...
                    psi_old(i, j+1) + psi_old(i, j-1) - h^2 * omega(i, j));
            end
        end
        error_psi = max(abs(psi(:) - psi_old(:)));
        % Update omega
        for i = 2:n-1
            for j = 2:m-1
                omega(i, j) = (1 - r_omega) * omega_old(i, j) + ...
                    r_omega * 0.25 * (omega_old(i+1, j) + omega_old(i-1, j) + ...
                    omega_old(i, j+1) + omega_old(i, j-1) - ...
                    (4 - h^2) * omega_old(i, j)) / (4 - h^2);
            end
        end
        error_omega = max(abs(omega(:) - omega_old(:)));
        error = max(error_psi, error_omega);
    end
    plot_Re10(psi);
end
%% Code snippet for plotting 'psi'
function plot_Re10(psi)
    % This function plots the streamfunction psi
    figure;
    contourf(psi', 20); % Transpose psi to get correct orientation
    colorbar;
    title('Streamfunction \psi for Re = 10');
    xlabel('x');
    ylabel('y');
end
댓글 수: 0
  Sam Chak
      
      
 2024년 8월 8일
        Hi @SANTHOSH
The code for plot_Re10() function is based on the video lecture of Prof. Jeffrey Chasnov. 
[psi, omega] = flow_around_cylinder_steady;
function [psi, omega] = flow_around_cylinder_steady
    Re=10; 
    %%%%% define the grid %%%%%
    n=101; m=101; % number of grid points
    N=n-1; M=m-1; % number of grid intervals
    h=pi/M; % grid spacing based on theta variable
    xi=(0:N)*h; theta=(0:M)*h; % xi and theta variables on the grid
    %%%%% Initialize the flow fields %%%%%
    psi=zeros(n,m);
    omega=zeros(n,m);
    psi(n,:)=exp(xi(n))*sin(theta(:)); % Write the free stream bc here
    %%%%% Set relax params, tol, extra variables %%%%%
    r_psi=1.8; % Set the relaxation parameter here, psi equation
    r_omega=0.9; % Set the relaxation parameter here, omega equation
    delta=1.e-08; % error tolerance
    error=2*delta; % initialize error variable
    %%%%% Add any additional variable definitions here %%%%%
    ...
    ...
    %%%%% Main SOR Loop %%%%%
    while (error > delta)
        psi_old = psi; omega_old = omega;
        for i=2:n-1
            for j=2:m-1
                psi(i,j)=(1-r_psi)*psi(i,j)+(r_psi/4)*(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+h*h*exp(2*xi(i))*omega(i,j));% Write psi equation here
            end
        end
        error_psi=max(abs(psi(:)-psi_old(:)));
        omega(1,:)=(psi(3,:) - 8*psi(2,:))/(2*h^2); % Write the boundary condition here
        for i=2:n-1
            for j=2:m-1
                omega(i,j)=(1-r_omega)*omega(i,j)+(r_omega/4)*(omega(i+1,j)+omega(i-1,j)+omega(i,j+1)+omega(i,j-1)+(Re/8)*((psi(i+1,j)-psi(i-1,j))*(omega(i,j+1)-omega(i,j-1))-(psi(i,j+1)-psi(i,j-1))*(omega(i+1,j)-omega(i-1,j)))); % Write omega equation here
            end
        end
        error_omega=max(abs(omega(:)-omega_old(:)));
        error=max(error_psi, error_omega);
    end
    %psi
    plot_Re10(psi);
    % figure(1)
    % contourf(xi,theta,psi.')
    % colorbar
    % figure(2)
    % contourf(xi,theta,omega.')
    % colorbar
end
%% Code by Prof. Jeffrey Chasnov
function plot_Re10(psi)
    % Reynolds number (design parameter)
    Re      = 10;
    % setup the xi-theta grid 
    n       = size(psi,1); 
    m       = size(psi,2);
    N       = n - 1; 
    M       = m - 1;
    h       = pi/M;     % grid spacing
    xi      = (0:N)*h; 
    theta   = (0:M)*h;
    [XI, THETA] = meshgrid(xi,theta);
    ... (Code is very long. Please refer to Prof. Jeffrey Chasnov's video lecture)
end
댓글 수: 6
  Cris LaPierre
    
      
 2024년 8월 9일
        Assuming this exercise is coming from Week 2 of Mathematics for Engineers: The Capstone Course, I can't duplicate the error. I'd suggest trying to submit your solution again.
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!











