Upload New File

Driven_Cavity_Solver.m

+function Driven_Cavity_Solver()
+    %% Driven Cavity, Finite Volume Method
+    % Paul Bataillie - paul.bataillie@student.isae-supaero.fr
+    % ---------------------------------------------------------------------
+    
+    close all
+    clc
+    clear
+    
+    %% Constants and parameters
+    Lx = 1;
+    Ly = 1;
+    rho = 1000;
+    mu = 1e-3;
+    
+    % Underrelaxation factors
+    e  = 0.5;
+    ep = 0.2;
+     
+    % Max iterations
+    numItr = 10000;
+    numItrP = 800;
+    dispResiduals = 20;
+    
+    % Number of cells 
+    Npx = 20;
+    
+    disp(append("INFO: Initializing Case ",num2str(Npx),"x",num2str(Npx),"."))
+    disp(append("INFO: Under-Relaxation Factors | e = ",num2str(e)," | ep = ",num2str(ep)))
+    
+    Npy = Npx;  % The pressure matrix (grid) is a square
+    
+    Nvx = Npx+2;  % The v-matrix contains two more columns and one more row
+    Nvy = Npx+1;  % than the pressure matrix (including dummy cells) 
+    
+    Nux = Npx+1;  % The u-matrix contains one more column and two more rows
+    Nuy = Npy+2;  % than the pressure matrix (including dummy cells) 
+    
+    % Cell size
+    dx = Lx/(Npx);
+    dy = Ly/(Npy);
+    
+    % Residual tolerances
+    RCE = 1e-5;
+    RME = 1e-7;
+    
+    %Reynolds number
+    Re = 1000;
+    Re1 = 1/Re;
+    
+    % Some variables
+    Dew = Re1*dy/dx;
+    Dns = Re1*dx/dy;
+    
+    % Initialization
+    u = zeros(Nuy,Nux);    
+    unew = zeros(Nuy,Nux); 
+    
+    v = zeros(Nvy,Nvx);
+    vnew = zeros(Nvy,Nvx);
+    
+    p  = zeros(Npy,Npx);
+    pp = zeros(Npy+2,Npx+2);
+    pp1 = zeros(Npy+2,Npx+2);
+    
+    auP = 0*ones(Nuy,Nux);
+    avP = 0*ones(Nvy,Nvx);
+    
+    uo = u;
+    vo = v;
+    
+    scrsz = get(0,'ScreenSize');
+    
+    %Boundray Conditions
+    unew(end,:) = 2 - unew(end-1,:);
+    u(end,:) = 2 - u(end-1,:);
+    
+    %Residuals Variables
+    continuity_res_scale_value = 0;
+    u_residual_vector = zeros(numItr,1);
+    v_residual_vector = zeros(numItr,1);
+    c_residual_vector = zeros(numItr,1);
+    
+    %% Start of the solution process
+    format long
+    
+    disp("INFO: Solving...")
+    
+    for i=1:numItr
+        %% Solve momentum equations
+        %% ------
+        % Mex
+        % Average velocities for the flux terms
+        uue = 0.5*(u(2:end-1,3:end) + u(2:end-1,2:end-1));
+        uuw = 0.5*(u(2:end-1,1:end-2) + u(2:end-1,2:end-1));
+        vun = 0.5*(v(2:end,2:end-2) + v(2:end,3:end-1));
+        vus = 0.5*(v(1:end-1,2:end-2) + v(1:end-1,3:end-1));
+        
+        %Flux terms
+        Fue = dy*uue;
+        Fuw = dy*uuw;
+        Fun = dx*vun;
+        Fus = dx*vus;
+        
+        % Coefficients
+        auE = Dew + max(0,-Fue);
+        auW = Dew + max(0,Fuw);
+        auN = Dns + max(0,-Fun);
+        auS = Dns + max(0,Fus);
+        auP = auE + auW + auN + auS + Fue + Fun - Fuw - Fus;
+    
+        % Solve MEx
+        unew(2:end-1,2:end-1) = (auE.*u(2:end-1,3:end) + auW.*u(2:end-1,1:end-2) + auN.*u(3:end,2:end-1) + auS.*u(1:end-2,2:end-1) - (p(1:end,2:end) - p(1:end,1:end-1))*dy) ./ (auP);
+        
+        % Update boundary conditions
+        unew(end,:) = 2 - unew(end-1,:);
+        unew(:,1) = 0;
+        unew(:,end) = 0;
+        unew(1,:) = - unew(2,:);
+        
+        %% ------
+        % Mey
+        % Average velocities for the flux terms
+        uve = 0.5*(u(2:end-2,2:end) + u(3:end-1,2:end));
+        uvw = 0.5*(u(2:end-2,1:end-1) + u(3:end-1,1:end-1));
+        vvn = 0.5*(v(2:end-1,2:end-1) + v(3:end,2:end-1));
+        vvs = 0.5*(v(1:end-2,2:end-1) + v(2:end-1,2:end-1));
+    
+        %Flux terms
+        Fve = dy*uve;
+        Fvw = dy*uvw;
+        Fvn = dx*vvn;
+        Fvs = dx*vvs;
+        
+        % Coefficients
+        avE = Dew + max(0,-Fve);
+        avW = Dew + max(0,Fvw);
+        avN = Dns + max(0,-Fvn);
+        avS = Dns + max(0,Fvs);
+        avP = avE + avW + avN + avS + Fve + Fvn - Fvw - Fvs;
+        
+        % Solve MEy
+        vnew(2:end-1,2:end-1) = (avE.*v(2:end-1,3:end) + avW.*v(2:end-1,1:end-2) + avN.*v(3:end,2:end-1) + avS.*v(1:end-2,2:end-1) - (p(2:end,1:end) - p(1:end-1,1:end))*dx) ./ (avP);
+        
+        % Update boundary conditions
+        vnew(1,:) = 0;
+        vnew(end,:) = 0;
+        vnew(:,1) = - vnew(:,2);
+        vnew(:,end) = - vnew(:,end-1);
+        
+        %% ------
+        % Pressure correction 
+        pp1 = zeros(Npy+2,Npx+2);
+    
+        apE = zeros(Npy,Npx); 
+        apN = zeros(Npy,Npx);
+        apW = zeros(Npy,Npx);
+        apS = zeros(Npy,Npx);
+        
+        apE(:,1:end-1) = dy^2./auP;
+        apE(:,end) = 0;
+    
+        apN(1:end-1,:) = dx^2./avP;
+        apN(end,:) = 0;
+    
+        apW(:,2:end) = dy^2./auP;
+        apW(:,1) = 0;
+    
+        apS(2:end,:) = dx^2./avP;
+        apS(1,:) = 0;
+    
+        ap = apE + apW + apN + apS;
+        bp = (-unew(2:end-1,2:end)+unew(2:end-1,1:end-1))*dy + (-vnew(2:end,2:end-1)+vnew(1:end-1,2:end-1))*dx;
+        
+        % Loop for p'
+        for j=1:numItrP
+            pp1(2:end-1,2:end-1) = ( apE.*pp(2:end-1,3:end) + apW.*pp(2:end-1,1:end-2) + apN.*pp(3:end,2:end-1) + apS.*pp(1:end-2,2:end-1) + bp ) ./ ap;
+            
+            % Underrelaxation
+            pp = ep*pp1 + (1-ep)*pp;
+    
+        end %...of p' loop  
+         
+        % Update p
+        p = p + pp(2:end-1,2:end-1);
+        
+        % Update u and v in internal points
+        u(2:end-1,2:end-1) = e * ( unew(2:end-1,2:end-1) + ( (dy) * (-pp(2:end-1,3:end-1) + pp(2:end-1,2:end-2)) )./(auP) ) + (1-e)*u(2:end-1,2:end-1);
+        v(2:end-1,2:end-1) = e * ( vnew(2:end-1,2:end-1) + ( (dx) * (-pp(3:end-1,2:end-1) + pp(2:end-2,2:end-1)) )./(avP) ) + (1-e)*v(2:end-1,2:end-1);
+       
+        % Update boundary conditions
+        u(end,:) = 2 - unew(end-1,:);
+        u(:,1) = 0;
+        u(:,end) = 0;
+        u(1,:) = - unew(2,:);
+        
+        v(1,:) = 0;
+        v(end,:) = 0;
+        v(:,1) = - vnew(:,2);
+        v(:,end) = - vnew(:,end-1);
+        
+        if (mod(i,dispResiduals)==0)||(i < 6)
+            %% Residuals
+            % U Momentum Sclaed Residual
+            num = abs(unew(2:end-1,2:end-1).*(auP) - ( (auE.*u(2:end-1,3:end) + auW.*u(2:end-1,1:end-2) + auN.*u(3:end,2:end-1) + auS.*u(1:end-2,2:end-1) - (p(1:end,2:end) - p(1:end,1:end-1))*dy) ));
+            den = abs(unew(2:end-1,2:end-1).*(auP));
+            r_u = sum(num,"all")/sum(den,"all");
+        
+            % V Momentum Scaled Residual
+            num = abs(vnew(2:end-1,2:end-1).*(avP) - ( (avE.*v(2:end-1,3:end) + avW.*v(2:end-1,1:end-2) + avN.*v(3:end,2:end-1) + avS.*v(1:end-2,2:end-1) - (p(2:end,1:end) - p(1:end-1,1:end))*dx) ));
+            den = abs(vnew(2:end-1,2:end-1).*(avP));
+            r_v = sum(num,"all")/sum(den,"all");
+    
+            % Continuity Scaled Residual (scaled as in Fluent, the scaled value is the largest absolute value of the continuity residual in the first five iterations)    
+            num = sum( abs( (u(2:end-1,2:end)-u(2:end-1,1:end-1))*dy + (v(2:end,2:end-1)-v(1:end-1,2:end-1))*dx ), "all");
+            if i < 6
+                continuity_res_scale_value = max(continuity_res_scale_value,num);
+            end
+            den = continuity_res_scale_value;        
+            r_c = num/den;
+            
+            % Monitor points
+            disp('----')
+            disp(append("numItr = ",num2str(i)," | U Residual = ",num2str(r_u)," | V Residual = ",num2str(r_v)," | C Residual = ",num2str(r_c)))
+            disp(append("Monitor Point | u = ",num2str(u(floor(Nux/4),floor(Nuy/4)))," | v = ",num2str(v(floor(Nux/4),floor(Nuy/4)))))
+    
+            u_residual_vector(i)=r_u;
+            v_residual_vector(i)=r_v;
+            c_residual_vector(i)=r_c;
+    
+            if (r_u<RME) && (r_v<RME) && (r_c<RCE)
+                break
+            end
+    
+        end    
+    
+    end
+    
+    %% Plotting 
+    
+    % Useful commands: meshgrid, streamslice, quiver, contour
+    
+    u_plot = zeros(Npy+2,Npx+2);
+    v_plot = zeros(Npy+2,Npx+2);
+    u_plot(2:end-1,2:end-1) = 0.5*(u(2:end-1,1:end-1) + u(2:end-1,2:end));
+    v_plot(2:end-1,2:end-1) = 0.5*(v(1:end-1,2:end-1) + v(2:end,2:end-1));
+    
+    u_plot(1,:) = 0;
+    u_plot(end,:) = 1;
+    u_plot(:,1) = 0;
+    u_plot(:,end) = 0;
+    
+    v_plot(1,:) = 0;
+    v_plot(end,:) = 0;
+    v_plot(:,1) = 0;
+    v_plot(:,end) = 0;
+    
+    x = zeros(1,Npx+2);
+    y = zeros(1,Npy+2);
+    
+    indx = (0:1/Npx:1) + 1/(Npx*2);
+    indy = (0:1/Npy:1) + 1/(Npy*2);
+    
+    x(1,2:end-1) = indx(1,1:end-1);
+    y(1,2:end-1) = indy(1,1:end-1);
+    x(1,end) = 1;
+    y(1,end) = 1;
+    
+    [X,Y] = meshgrid(x,y);
+    
+    figure
+    streamslice(X,Y,u_plot,v_plot,5);
+    title('Velocity Streamlines')
+    ylabel('$\frac{y}{L}$','interpreter','latex')
+    xlabel('$\frac{x}{L}$','interpreter','latex')
+    set(gca,'ylim',[0 1]);
+    axis tight
+    
+    u_y_line = u_plot(:,ceil((Npx+2)/2));
+    v_x_line = v_plot(ceil((Npx+2)/2),:);
+    
+    format short
+end