Matlab code for pipe flow
Hi All,
I have written a MATLAB code for a 2D lid driven cavity problem, and it works fast and well, the results are consistent with experimental data. But this code does not work for the duct flow, which is also a classic example in many references. I wonder it is due to the change of the definition of boundary conditions or the scheme itself. What is the main issue behind this. Thanks. 
Well done writing your own Matlab code!
If all is well (you are correctly solving NavierStokes equations) changing Lid Driven Cavity flow is as 'easy' as changing the Boundary Conditions. Do you make use of a staggered (velocity defined on cell faces, pressure in cell center) grid or a collocated arrangement (velocity and pressure defined in cell center)? You should clearly specify the BCs you want to use for your flow. The following seems logical: Inflow: uniform inflow profile U = 1, dvdx=0 (freeslip) and dpdx=0 (neumann) Walls: U = 0, V = 0, dpdy=0 Ouflow: free outflow dudx = 0, dvdx = 0 and dpdx = 0 You could also use dudx = dvdy at the outflow (from continuity equation). You can also try freeslip BCs for the velocity at the walls, which will better represent infinite domain then the velocity=0. Also setting one of the pressures to a reference value could be usefull (p = 0 at outflow for example). The way these BCs translate to your code depend on whether the grid is staggered or collocated. Also make sure you correctly alter the matrices you use for your computation besides altering the boundary values. Discretization schemes should not change when you differ the BCs. However the matrices used for multiplication change depending on the boundary conditions, the value at matrix location (1,1) and (n,n). The value should be equal to 1 for dudx equal to zero and 2 for u equal to zero (check this for all the BCs, this example is for 2nd differences in space). Well good luck and enjoy, hopefully you can use my comments. Note: you should interpret the third dimension of a 2D simulation as infinite, so instead of pipe flow you are simulation flow through 2 horizontal infinite plates (in that third dimension). 
Thanks, Vincent.
I am using the fractional step method (predictor and corrector) on a staggered grid. The boundary conditions I defined for the pipe flow are as follows: Inlet: uniform u velocity, u = 1 (nondimensionalized), v = 0, and dpdx = 0 Walls: Walls: u = 0, v= 0, dpdy=0 Outlet: dudx = 0, dvdx = 0, and p = 0 The scheme is fully explicit, so no additional intermediate boundary conditions are necessary. However, when it comes to the Pressure Poisson Equation, the computation process slows down a lot and gives impractical results. Also, the location (1,1) and (n,n) never enters the computation in my code, since I use a staggered grid. As to the note, what I really wanna simulate is the cylindrical pipe flow. But your comment is really important, since it determines the problem to be a 2D or axisymmetric. I am still tracking the problem, it is really time consuming. Thanks again, and if you have further idea about what is going on with this crazy problem, please let me know. 
Hey there,
your BCs look just fine, just be care of setting the boundary conditions. On a staggered grid setting u = 0 (wall condition in the ydirection) is done as setting u(bc) = u(end) so your boundary condition will the 1*the value of u close to the wall, due to the fact your u cell is not defined on the wall. Setting the boundary condition dudx can be done by setting u(bc) = u(end) dx*p(end), or u(bc)=u(end) since p(end) = 0. However you should do this step after you pressure correction step (or befor you determine the prediction velocity field). Let me clearfy a point I tried to make in my post but that didn't come across: "Discretization schemes should not change when you differ the BCs. However the matrices used for multiplication change depending on the boundary conditions, the value at matrix location (1,1) and (n,n). The value should be equal to 1 for dudx equal to zero and 2 for u equal to zero (check this for all the BCs, this example is for 2nd differences in space)." In Matlab you should determine your velocity/pressure fields by matrix multiplications. For instance taking the centered second derivative of V can be done by: Creating a 2nd difference matrix U: U = 2*eye(n)  diag(ones(n1,1),1)  diag(ones(n1,1),1) (2 on the diagonal and 1 on the two off diagonals). And multiplying this with your velocity matrix V: d^2V/dx^2 = 1/(dx^2) U*V But we should be very aware of the first and last line of this matrix U. Namely using this form the first element is determined by 2*V(i,j)  1*V(i,j+1), which is fine for dirichlet conditions on a point (u in x direction). But should to 1 (neumann) or 3 (dirichlet between 2 points) in order to function correctly. 2D flow of 2 parallel plates is not the same as flow through a pipe, take for instance permeability of the pipe/plates. In case of a pipe K = R^2/8 where for 2 plates K = 'R'^2/3, so be carefull interpreting your results from the 2D simulation. I would also suggest you go to cylindrical coordinates if you really want to simulate pipeflow. Good luck! 
Thanks for clear explanation.
Just for curiosity's sake, have you tried to implement this problem? Or you give suggestions from your experience of CFD work? :) 
I'm doing an internship at a company at the moment and was asked to program in Matlab. I had previous experience with Fortran but had to relearn some work so I started on 2D DNS of flow between 2 plates.
Now I'm working on incorperating the Volume Averaged Navier Stokes (VANS) equations and using an Immersed Boundary Method (IBM) to introduce my geometry. So far I'm really enjoying Matlab, I can get things done at lot quicker (and with less coding) than in Fortran (at the cost of speed/efficiency). One program to do all your work (figures etc) is very nice as well. Regards, Vincent 
Hi, I have found a MATLAB code online for the liddriven cavity problem, and trying to modify it to see if this code will work for a channel flow or Couette flow. Do you mind having a look at this one? I think this will help revise my own code as well.
Best 
Sure, I will gladly help you out. Just post it on the forum so that other people that are interested can follow the discussion.
Regards, Vincent 
Here it is. there are 2 mfiles.
1. pcolm.m (it is splitted in to 2 parts, since too long) %%% pcolm  a simplified version of pcol.f code (Peric, 1997) %%% %%% This version is writen for matlab. %%% %%% Gabriel Usera  Feb/2007 %%% %%% This version is setup for LidDriven Cavity flow. %%% %%% To set up other flows B.C. require modification. %%% %%% Typical user input is covered in lines 43 to 90, %%% including all numerical parameters, fluid properties %%% and grid specification (lines 89 and 90). %%% %%% Boundary Conditions are specified in lines : %%% 126  Initial B.C. setup for Velocity %%% 288332  Complete B.C. for Velocity %%% 431435  B.C. for mass balance/pressure equation %%% %%% Please keep this version 'as it is' and rename %%% the versions that you modify ('poclm9538.m' or %%% something like that). %%% %%% To run : %%% >> [X,Y,XC,YC,FX,FY,U,V,P,F1,F2]=pcolm; %%% %%% To plot : %%% >> plotm; %%% %%% function [X,Y,XC,YC,FX,FY,U,V,P,F1,F2]=pcolm; global LTIME LTEST % Logic variables, Steady/Transient and Testing global IPR JPR IU IV IP % (IPR,JPR) Pressure reference. IU,IV,IP Equation tags. global NI NJ NIM NJM NIJ LI % Domain dimensions and indexing vector LI=(I1)*NJ global X Y XC YC FX FY % Domain coordinates and interpolation factors global DENSIT VISC GDS DT DTR SMALL ALFA ULID % Fluid properties, Blending factor, Time step,... global U V P PP F1 F2 DPX DPY U0 V0 % Variables: velocities, pressure, mass flux, pressure gradient... global AE AW AS AN AP APR SU SV % Coefficient matrices global SOR URF NSW RESOR % Iteration control parameters %%% INPUT DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% ITIM=0; % Set iteration counter to 0 TIME=0.; % Set initial time to 0. % LTEST=0; % Logic Flag : Testing ? LTEST=1, otherwise LTEST=0. LTIME=0; % Logic Flag : Stationary (0) or Transient (1) ? % MAXIT=500; % Maximum number of outer iterations in each time step IMON=5; % Index (I) for monitoring location JMON=5; % Index (J) for monitoring location IPR=2; % Index (I) for pressure reference point JPR=2; % Index (J) for pressure reference point SORMAX=1e4; % Residual level for stopping outer iterations SLARGE=1e+3; % Residual level for divergence ALFA=0.92; % Parameter for linear solver SIPSOL % DENSIT=1.0; % Fluid density VISC= 1.e2; % Fluid viscosity ULID= 1.; % Lid velcity for lid driven cavity % UIN=0.; % Initial value for U velocity (usually 0.) VIN=0.; % Initial value for V velocity (usually 0.) PIN=0.; % Initial valur for Pressure (usually 0.) % ITST=1; % Number of time steps (1 for steady flow, LTIME=0) DT=1.00; % Time step in seconds (meaningless for steady flow) % IU=1; % Tag for U equation IV=2; % Tag for V equation IP=3; % Tag for P equation % URF(IU)=0.8; % Under relaxation parameter for U URF(IV)=0.8; % Under relaxation parameter for V URF(IP)=0.2; % Under relaxation parameter for P % SOR(IU)=0.2; % Ratio of residual reduction for linear solver : U SOR(IV)=0.2; % Ratio of residual reduction for linear solver : V SOR(IP)=0.2; % Ratio of residual reduction for linear solver : P % NSW(IU)=1; % Maximum number of linear solver iterations : U NSW(IV)=1; % Maximum number of linear solver iterations : V NSW(IP)=6; % Maximum number of linear solver iterations : P % GDS=1.0; % UDS  CDS Blending for U,V % %%% GRID DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% X=[0:.025:1,1]'; % X  cell volume coordinates (only tested uniform spacing) Y=[0:.025:1,1]'; % Y  cell volume coordinates (only tested uniform spacing) % NI=length(X); % Number of points in X ( + 1 ) NJ=length(Y); % Number of points in Y ( + 1 ) NIM=NI1; % NI  1 NJM=NJ1; % NJ  1 NIJ=NI*NJ; % Total number of points %%% Indexing vector for k=(i1)*NJ+j indexing style LI=([1:NI]1)*NJ; % %%% X cell center coordinate XC=X; XC(2:NIM)=0.5*(X(1:NIM1)+X(2:NIM)); %%% Y cell center coordinate YC=Y; YC(2:NJM)=0.5*(Y(1:NJM1)+Y(2:NJM)); %%% Interpolation factors (in uniform grid FX=FY=0.5 except at boundaries) FX(1:NIM)=(X(1:NIM)XC(1:NIM))./(XC(2:NI)XC(1:NIM)); FX(NI)=0.; FY(1:NJM)=(Y(1:NJM)YC(1:NJM))./(YC(2:NJ)YC(1:NJM)); FY(NJ)=0.; %%% SET SOME CONTROL VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SMALL=1.e15; % Small number GREAT=1.e+15; % Big number DTR=1./DT; % Reciprocal of time step RESOR=zeros(3,1); % Initialize array to store residuals %%% Monitoring location IJMON=LI(IMON)+JMON; %%% INITIALIZE FIELDS and ARRAYS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [U,V,P,PP,U0,V0,F1,F2,DPX,DPY]=deal(zeros(NIJ,1)); [AP,AE,AW,AN,AS,APR,SU,SV]=deal(zeros(NIJ,1)); %%% FIXED BOUNDARY AND INITIAL CONDITIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% North wall velocity for I=2:NIM, U(LI(I)+NJ)=ULID; end; %%% Intial conditions for I=2:NIM, for IJ=LI(I)+2:LI(I)+NJM, U(IJ)=UIN; V(IJ)=VIN; P(IJ)=PIN; U0(IJ)=UIN; V0(IJ)=VIN; end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%% TIME LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% ITIMS=ITIM+1; % Update iteration counter to first iteration ITIME=ITIM+ITST; % Final number of iterations % for ITIM=ITIMS:ITIME, % Set time loop TIME=TIME+DT; % Update time % if LTIME, U0=U; V0=V; end; % Shift solutions in time % %%% HEADDING FOR THIS TIME STEP disp([' ']); disp(['************************************************* *************']); disp(['TIME=',num2str(TIME,'%0.2E%')]); disp([' ']); disp(['IT.RES(U)RES(V)RES(P)UMONVMONPMON']); % % %%% OUTER ITERATIONS LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ITER=1:MAXIT CALCUV; % Call CALCUV routine to update U and V CALCP; % Call CALCP routine to update P,F1,F2 and U,V again % % Display information about residuals and monitor point values disp([num2str(ITER),' ',num2str(RESOR(IU),'%0.2E'),' ',num2str(RESOR(IV),'%0.2e'),... ' ',num2str(RESOR(IP),'%0.2e'),' ',num2str(U(IJMON) ,'%0.2e'),... ' ',num2str(V(IJMON) ,'%0.2e'),' ',num2str(P(IJMON) ,'%0.2e')]); % % Check convergence SOURCE=max(RESOR([IU,IV,IP])); if SOURCE>SLARGE, break; end; if SOURCE<SORMAX, break; end; end; % % Return/Break if DIVERGING if SOURCE>SLARGE, disp('DIVERGING'); return; end; end; % disp([' ']); disp(['CALCULATION FINISHED  SEE RESULTS IN [X,Y,XC,YC,FX,FY,U,V,P,F1,F2]']); % return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% function CALCUV global LTIME LTEST global IPR JPR IU IV IP global NI NJ NIM NJM NIJ LI global X Y XC YC FX FY global DENSIT VISC GDS DT DTR SMALL ALFA ULID global U V P PP F1 F2 DPX DPY U0 V0 global AE AW AS AN AP APR SU SV global SOR URF NSW RESOR %%% INITIALIZE COEFFICIENTS AND SOURCE TERMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [AP,AE,AW,AN,AS,SU,SV,APR]=deal(zeros(NIJ,1)); %%% FLUXES THROUGH INTERNAL EAST CV FACES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM1, FXE=FX(I); % Interpolation Factor FXP=1.FXE; % DXPE=XC(I+1)XC(I); % Distance P>E % for J=2:NJM, IJ=LI(I)+J; % Index IJ to P IJE=IJ+NJ; % Index IJE to E % S=Y(J)Y(J1); % Cell face 'area' D=VISC*S/DXPE; % Coefficient for diffusive flux % CE=min(F1(IJ),0.); % Mass fluxes, upwind from E CP=max(F1(IJ),0.); % Mass fluxes, upwind from P % FUUDS=CP*U(IJ)+CE*U(IJE); % Explicit convective fluxes UDS, U FVUDS=CP*V(IJ)+CE*V(IJE); % Explicit convective fluxes UDS, V FUCDS=F1(IJ)*(U(IJE)*FXE+U(IJ)*FXP); % Explicit convective fluxes CDS, U FVCDS=F1(IJ)*(V(IJE)*FXE+V(IJ)*FXP); % Explicit convective fluxes CDS, V % AE(IJ )=+CED; % Coefficient for E in P AW(IJE)=CPD; % Coefficient for W in E % SU(IJ )=SU(IJ )+GDS*(FUUDSFUCDS); % Source term for P, U SU(IJE)=SU(IJE)GDS*(FUUDSFUCDS); % Source term for E, U SV(IJ )=SV(IJ )+GDS*(FVUDSFVCDS); % Source term for P, V SV(IJE)=SV(IJE)GDS*(FVUDSFVCDS); % Source term for E, V end; end; %%% FLUXES THROUGH INTERNAL NORTH CV FACES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for J=2:NJM1, FYN=FY(J); % Interpolation Factor FYP=1.FYN; % DYPN=YC(J+1)YC(J); % Distance P>N % for I=2:NIM, IJ=LI(I)+J; % Index IJ to P IJN=IJ+1; % Index IJN to N % S=X(I)X(I1); % Cell face 'area' D=VISC*S/DYPN; % Coefficient for diffusive flux % CN=min(F2(IJ),0.); % Mass fluxes, upwind from E CP=max(F2(IJ),0.); % Mass fluxes, upwind from P % FUUDS=CP*U(IJ)+CN*U(IJN); % Explicit convective fluxes UDS, U FVUDS=CP*V(IJ)+CN*V(IJN); % Explicit convective fluxes UDS, V FUCDS=F2(IJ)*(U(IJN)*FYN+U(IJ)*FYP); % Explicit convective fluxes CDS, U FVCDS=F2(IJ)*(V(IJN)*FYN+V(IJ)*FYP); % Explicit convective fluxes CDS, V % AN(IJ )=+CND; % Coefficient for E AS(IJN)=CPD; % Coefficient for W % SU(IJ )=SU(IJ )+GDS*(FUUDSFUCDS); % Source term for P, U SU(IJN)=SU(IJN)GDS*(FUUDSFUCDS); % Source term for N, U SV(IJ )=SV(IJ )+GDS*(FVUDSFVCDS); % Source term for P, V SV(IJN)=SV(IJN)GDS*(FVUDSFVCDS); % Source term for N, V end; end; %%% VOLUME INTEGRALS (SOURCE TERMS) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM, % Loop I direction DX=X(I)X(I1); % Cell length for J=2:NJM, % Loop J direction DY=Y(J)Y(J1); % Cell heigth VOL=DX*DY; % Cell Volume IJ=LI(I)+J; % IJ index to P % PE=P(IJ+NJ)*FX(I )+P(IJ )*(1.FX(I )); % Pressure at 'e' PW=P(IJ )*FX(I1)+P(IJNJ)*(1.FX(I1)); % Pressure at 'w' PN=P(IJ+ 1)*FY(J )+P(IJ )*(1.FY(J )); % Pressure at 'n' PS=P(IJ )*FY(J1)+P(IJ 1)*(1.FY(J1)); % Pressure at 's' DPX(IJ)=(PEPW)/DX; % Pressure Xgradient DPY(IJ)=(PNPS)/DY; % Pressure Ygradient SU(IJ)=SU(IJ)DPX(IJ)*VOL; % Pressure term  U SV(IJ)=SV(IJ)DPY(IJ)*VOL; % Pressure term  V % if LTIME, % Unsteady ? APT=DENSIT*VOL*DTR; % Coefficient SU(IJ)=SU(IJ)+APT*U0(IJ); % Source term  U SV(IJ)=SV(IJ)+APT*V0(IJ); % Source term  V AP(IJ)=AP(IJ)+APT; % AP Coef. end; end; end; %%% BOUNDARY CONDITIONS U,V %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SOUTH BOUNDARY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% for I=2:NIM, % Loop I direction IJ=LI(I)+2; IJB=IJ1; % IJ,IJB to South Boundary (J=2) U(IJB)=0.; % Set South wall U velocity V(IJB)=0.; % Set South wall V velocity D=VISC*(X(I)X(I1))/(YC(2)YC(1)); % Diffusion coef. AP(IJ)=AP(IJ)+D; % Coef for P SU(IJ)=SU(IJ)+D*U(IJB); % Sourec term  U SV(IJ)=SV(IJ)+D*V(IJB); % Sourec term  V end; %%% NORTH BOUNDARY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% for I=2:NIM, IJ=LI(I)+NJM; IJB=IJ+1; U(IJB)=ULID; V(IJB)=0.; D=VISC*(X(I)X(I1))/(YC(NJ)YC(NJM)); AP(IJ)=AP(IJ)+D; SU(IJ)=SU(IJ)+D*U(IJB); SV(IJ)=SV(IJ)+D*V(IJB); end; %%% WEST BOUNDARY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% for J=2:NJM, IJ=LI(2)+J; IJB=IJNJ; U(IJB)=0.; V(IJB)=0.; D=VISC*(Y(J)Y(J1))/(XC(2)XC(1)); AP(IJ)=AP(IJ)+D; SU(IJ)=SU(IJ)+D*U(IJB); SV(IJ)=SV(IJ)+D*V(IJB); end; %%% EAST BOUNDARY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% for J=2:NJM, IJ=LI(NIM)+J; IJB=IJ+NJ; U(IJB)=0.; V(IJB)=0.; D=VISC*(Y(J)Y(J1))/(XC(NI)XC(NIM)); AP(IJ)=AP(IJ)+D; SU(IJ)=SU(IJ)+D*U(IJB); SV(IJ)=SV(IJ)+D*V(IJB); end; %%% SOLVE EQUATIONS FOR U AND V %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%% UNDER RELAXATION, SOLVING FOR UVELOCITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM, for IJ=LI(I)+2:LI(I)+NJM, APR(IJ)=AP(IJ); AP(IJ)=(APR(IJ)AE(IJ)AW(IJ)AN(IJ)AS(IJ))/URF(IU); SU(IJ)=SU(IJ)+(1.URF(IU))*AP(IJ)*U(IJ); end; end; %%% SOLVE for U U=SIPSOL(U,IU); %%% UNDER RELAXATION, SOLVING FOR VVELOCITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM, for IJ=LI(I)+2:LI(I)+NJM, AP(IJ)=(APR(IJ)AE(IJ)AW(IJ)AN(IJ)AS(IJ))/URF(IV); SU(IJ)=SV(IJ)+(1.URF(IV))*AP(IJ)*V(IJ); APR(IJ)=1./AP(IJ); end; end; %%% SOLVE for V V=SIPSOL(V,IV); return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% 
function CALCP;
global LTIME LTEST global IPR JPR IU IV IP global NI NJ NIM NJM NIJ LI global X Y XC YC FX FY global DENSIT VISC GDS DT DTR SMALL ALFA ULID global U V P PP F1 F2 DPX DPY U0 V0 global AE AW AS AN AP APR SU SV global SOR URF NSW RESOR %%% INITIALIZE COEFFICIENTS AND SOURCE TERMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [AP,AE,AW,AN,AS,SU]=deal(zeros(NIJ,1)); %%% EAST CV FACES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% for I=2:NIM1, % Loop I direction DXPE=XC(I+1)XC(I); % Distance P>E FXE=FX(I); % Interpolation factor FXP=1.FXE; % Interpolation factor % for J=2:NJM, % Loop J direction IJ=LI(I)+J; % Index IJ to P IJE=IJ+NJ; % Index IJE to E S=Y(J)Y(J1); % Cell face 'area' VOLE=DXPE*S; % Volume between P and E D=DENSIT*S; % Coefficient % DPXEL=0.5*(DPX(IJE)+DPX(IJ)); % Interpolated gradient in 'e' UEL=U(IJE)*FXE+U(IJ)*FXP; % Interpolated Uvelocity in 'e' APUE=APR(IJE)*FXE+APR(IJ)*FXP; % Interpolated 1/AP coeff. in 'e' % DPXE=(P(IJE)P(IJ))/DXPE; % Gradient in 'e', compact aprox. UE=UELAPUE*VOLE*(DPXEDPXEL); % Corrected U velocity in 'e' F1(IJ)=D*UE; % Mass flux through 'e' face % AE(IJ)=D*APUE*S; % Coefficient for E in P AW(IJE)=AE(IJ); % Coefficient for W in E end; end; %%% NORTH CV FACES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% for J=2:NJM1, % Loop J direction DYPN=YC(J+1)YC(J); % Distance P>N FYN=FY(J); % Interpolation factor FYP=1.FYN; % Interpolation factor % for I=2:NIM, % Loop I direction IJ=LI(I)+J; % Index IJ to P IJN=IJ+1; % Index IJN to N S=X(I)X(I1); % Cell face 'area' VOLN=DYPN*S; % Volume between P and N D=DENSIT*S; % Coefficient % DPYNL=0.5*(DPY(IJN)+DPY(IJ)); % Interpolated gradient in 'n' VNL=V(IJN)*FYN+V(IJ)*FYP; % Interpolated Vvelocity in 'n' APVN=APR(IJN)*FYN+APR(IJ)*FYP; % Interpolated 1/AP coeff. in 'n' % DPYN=(P(IJN)P(IJ))/DYPN; % Gradient in 'n', compact aprox. VN=VNLAPVN*VOLN*(DPYNDPYNL); % Corrected V velocity in 'n' F2(IJ)=D*VN; % Mass flux through 'n' face % AN(IJ)=D*APVN*S; % Coefficient for N in P AS(IJN)=AN(IJ); % Coefficient for S in N end; end; %%% SINCE ALL BOUNDARIES ARE ZERO MASS FLUX BOUNDARIES (WALLS), %%%%% %%% WE HAVE NEWMANN BOUNDARY CONDITIONS FOR PRESSURE CORRECTION %%%%% %%% (NULL GRADIENT). NO SPECIAL TREATMEN IS REQUIRED. %%%%% %%% FOR THE CASE OF INLETS AND OUTLETS, MASS FLUXES AT THE BOUNDARIES %%%%% %%% NEED TO BE COMPUTED HERE. %%%%% %%% SOURCE TERM AND COEFFICIENT FOR NODE P %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUM=0.; % Initialize SUM for I=2:NIM, % Loop I direction for IJ=LI(I)+2:LI(I)+NJM, % Loop J direction SU(IJ)=F1(IJNJ)F1(IJ)+F2(IJ1)F2(IJ); % Flux imbalance AP(IJ)=(AE(IJ)+AW(IJ)+AN(IJ)+AS(IJ)); % Coefficient for P SUM=SUM+SU(IJ); % Global flux imbalance PP(IJ)=0.; % Initialize PP end; end; if LTEST, disp(['SUM=',num2str(SUM)]); end; % If testing display SUM %%% solve for PP PP=SIPSOL(PP,IP); %%% EXTRAPOLATE BOUNDARY VALUES FOR PP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SOUTH AND NORTH BOUNDARIES for I=2:NIM, IJ=LI(I)+1; PP(IJ)=PP(IJ+1)+(PP(IJ+1)PP(IJ+2))*FY(2); IJ=LI(I)+NJ; PP(IJ)=PP(IJ1)+(PP(IJ1)PP(IJ2))*(1.FY(NJM1)); end; %%% WEST AND EAST BOUNDARIES NJ2=2*NJ; for J=2:NJM, IJ=LI(1)+J; PP(IJ)=PP(IJ+NJ)+(PP(IJ+NJ)PP(IJ+NJ2))*FX(2); IJ=LI(NI)+J; PP(IJ)=PP(IJNJ)+(PP(IJNJ)PP(IJNJ2))*(1.FX(NIM1)); end; %%% REFERENCE PRESSURE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% IJPREF=LI(IPR)+JPR; PP0=PP(IJPREF); %%% CORRECT EAST MASS FLUXES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM1 for IJ=LI(I)+2:LI(I)+NJM, F1(IJ)=F1(IJ)+AE(IJ)*(PP(IJ+NJ)PP(IJ)); end; end; %%% CORRECT NORTH MASS FLUXES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM for IJ=LI(I)+2:LI(I)+NJM1, F2(IJ)=F2(IJ)+AN(IJ)*(PP(IJ+1)PP(IJ)); end; end; %%% CORRECT PRESSURE AND VELOCITIES AT CELL CENTER %%%%%%%%%%%%%%%%%%%%%%%% for I=2:NIM, DX=X(I)X(I1); % for J=2:NJM, IJ=LI(I)+J; DY=Y(J)Y(J1); % PPE=PP(IJ+NJ)*FX(I )+PP(IJ )*(1.FX(I )); % Pressure at 'e' PPW=PP(IJ )*FX(I1)+PP(IJNJ)*(1.FX(I1)); % Pressure at 'w' PPN=PP(IJ+1 )*FY(J )+PP(IJ )*(1.FY(J )); % Pressure at 'n' PPS=PP(IJ )*FY(J1)+PP(IJ1 )*(1.FY(J1)); % Pressure at 's' % U(IJ)=U(IJ)(PPEPPW)*DY*APR(IJ); % Uvelocity corrected V(IJ)=V(IJ)(PPNPPS)*DX*APR(IJ); % Vvelocity corrected P(IJ)=P(IJ)+URF(IP)*(PP(IJ)PP0); % Pressure corrected end; end; %%% EXTRAPOLATE BOUNDARY VALUES FOR P %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SOUTH AND NORTH BOUNDARIES for I=2:NIM, IJ=LI(I)+1; P(IJ)=P(IJ+1)+(P(IJ+1)P(IJ+2))*FY(2); IJ=LI(I)+NJ; P(IJ)=P(IJ1)+(P(IJ1)P(IJ2))*(1.FY(NJM1)); end; %%% WEST AND EAST BOUNDARIES NJ2=2*NJ; for J=2:NJM, IJ=LI(1)+J; P(IJ)=P(IJ+NJ)+(P(IJ+NJ)P(IJ+NJ2))*FX(2); IJ=LI(NI)+J; P(IJ)=P(IJNJ)+(P(IJNJ)P(IJNJ2))*(1.FX(NIM1)); end; return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% function FI=SIPSOL(FI,IFI); global LTEST global NI NJ NIM NJM NIJ LI global SMALL ALFA global AE AW AS AN AP SU global SOR NSW RESOR %%% INITIALIZE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% [UE,UN,RES,LW,LS,LPR]=deal(zeros(NIJ,1)); %%% COEFFICIENTS OF UPPER AND LOWER TIRANGULAR MATRICES %%%%%%%%%%%%%%%%%%% for I=2:NIM, for IJ=LI(I)+2:LI(I)+NJM, LW(IJ)=AW(IJ)/(1.+ALFA*UN(IJNJ)); LS(IJ)=AS(IJ)/(1.+ALFA*UE(IJ1)); P1=ALFA*LW(IJ)*UN(IJNJ); P2=ALFA*LS(IJ)*UE(IJ1); LPR(IJ)=1./(AP(IJ)+P1+P2LW(IJ)*UE(IJNJ)LS(IJ)*UN(IJ1)); UN(IJ)=(AN(IJ)P1)*LPR(IJ); UE(IJ)=(AE(IJ)P2)*LPR(IJ); end; end; %%% INNER ITERATIONS LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for L=1:NSW(IFI), RESL=0.; % %%% CALCULATE RESIDUAL AND OVERWRITE IT BY INTERMEDIATE VECTOR%%%%%%%%% for I=2:NIM, for IJ=LI(I)+2:LI(I)+NJM, RES(IJ)=SU(IJ)AN(IJ)*FI(IJ+1)AS(IJ)*FI(IJ1)AE(IJ)*FI(IJ+NJ)... AW(IJ)*FI(IJNJ)AP(IJ)*FI(IJ); RESL=RESL+abs(RES(IJ)); RES(IJ)=(RES(IJ)LS(IJ)*RES(IJ1)LW(IJ)*RES(IJNJ))*LPR(IJ); end; end; % %%% STORE INITIAL RESIDUAL SUM FOR CHECKING CONVERGENCE OF OUTER IT. %% if L==1, RESOR(IFI)=RESL; end; RSM=RESL/(RESOR(IFI)+SMALL); % %%% BACK SUBSTITUTION AND CORRECTION for I=NIM:1:2, for IJ=LI(I)+NJM:1:LI(I)+2, RES(IJ)=RES(IJ)UN(IJ)*RES(IJ+1)UE(IJ)*RES(IJ+NJ); FI(IJ)=FI(IJ)+RES(IJ); end; end; % %%% CHECK CONVERGENCE if LTEST, disp([num2str(L),' INNER ITER., RESL=',num2str(RESL)]); end; if RSM<SOR(IFI), return; end; end; % return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% 
2. plotm.m
%%% Run pcolm [X,Y,XC,YC,FX,FY,U,V,P,F1,F2]=pcolm; %%% Get size NI=length(X); NJ=length(Y); sz=[NJ,NI]; %%% Quiver velocity field figure; quiver(XC,YC,reshape(U,sz),reshape(V,sz)); axis image; %%% Streamlines figure; h=streamline(XC',YC,reshape(U,sz),reshape(V,sz),ra nd(50,1),rand(50,1),[0.1,1000]); axis image; axis([0,1,0,1]); 
Wow, thats alot of code. One thing you should not do in matlab is use for loops, it's terrribly slow and using matrix multiplication will result in a much cleaner code. So I'm going to make a counter proposition on the code to use, it's alot more clear and it has a good document which explains the code.
PS: I take no credit whatsoever for this code, it's written by Benjamin Seibold (MIT). The document can be found at: http://wwwmath.mit.edu/cse/codes/mi...vierstokes.pdf The m file can be found at: http://wwwmath.mit.edu/cse/codes/mi...navierstokes.m I posted the code here in order to discuss the code and on how we can modify it (it's Lid Cavity atm). I suggest you first read the document and code a couple times to fully understand what's going on. Take care! Vincent function mit18086_navierstokes %MIT18086_NAVIERSTOKES % Solves the incompressible NavierStokes equations in a % rectangular domain with prescribed velocities along the % boundary. The solution method is finite differencing on % a staggered grid with implicit diffusion and a Chorin % projection method for the pressure. % Visualization is done by a colormapisoline plot for % pressure and normalized quiver and streamline plot for % the velocity field. % The standard setup solves a lid driven cavity problem. % 07/2007 by Benjamin Seibold % http://wwwmath.mit.edu/~seibold/ % Feel free to modify for teaching and learning. % Re = 1e2; % Reynolds number dt = 1e2; % time step tf = 4e0; % final time lx = 1; % width of box ly = 1; % height of box nx = 90; % number of xgridpoints ny = 90; % number of ygridpoints nsteps = 10; % number of steps with graphic output % nt = ceil(tf/dt); dt = tf/nt; x = linspace(0,lx,nx+1); hx = lx/nx; y = linspace(0,ly,ny+1); hy = ly/ny; [X,Y] = meshgrid(y,x); % % initial conditions U = zeros(nx1,ny); V = zeros(nx,ny1); % boundary conditions uN = x*0+1; vN = avg(x)*0; uS = x*0; vS = avg(x)*0; uW = avg(y)*0; vW = y*0; uE = avg(y)*0; vE = y*0; % Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hx^2+... [uW;zeros(nx3,ny);uE]/hy^2); Vbc = dt/Re*([vS' zeros(nx,ny3) vN']/hx^2+... [2*vW(2:end1);zeros(nx2,ny1);2*vE(2:end1)]/hy^2); fprintf('initialization') Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye (nx)); Lp(1,1) = 3/2*Lp(1,1); perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp'; Lu = speye((nx1)*ny)+dt/Re*(kron(speye(ny),K1(nx1,hx,2))+... kron(K1(ny,hy,3),speye(nx1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru'; Lv = speye(nx*(ny1))+dt/Re*(kron(speye(ny1),K1(nx,hx,3))+... kron(K1(ny1,hy,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; Lq = kron(speye(ny1),K1(nx1,hx,2))+kron(K1(ny1,hy,2),speye(nx1)); perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; fprintf(', time loop\n20%%40%%60%%80%%100%%\n') for k = 1:nt % treat nonlinear terms gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ue = [2*uS'Ue(:,1) Ue 2*uN'Ue(:,end)]; Ve = [vS' V vN']; Ve = [2*vWVe(1,:);Ve;2*vEVe(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Vagamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Vagamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end1)); Ud = diff(Ue(:,2:end1))/2; Va = avg(Ve(2:end1,:)')'; Vd = diff(Ve(2:end1,:)')'/2; U2x = diff(Ua.^2gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2gamma*abs(Va).*Vd)')'/hy; U = Udt*(UVy(2:end1,:)+U2x); V = Vdt*(UVx(:,2:end1)+V2y); % implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny1); % pressure correction rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); p(perp) = Rp\(Rpt\rhs(perp)); P = reshape(p,nx,ny); U = Udiff(P)/hx; V = Vdiff(P')'/hy; % visualization if floor(25*k/nt)>floor(25*(k1)/nt), fprintf('.'), end if k==1floor(nsteps*k/nt)>floor(nsteps*(k1)/nt) % stream function rhs = reshape(diff(U')'/hydiff(V)/hx,[],1); q(perq) = Rq\(Rqt\rhs(perq)); Q = zeros(nx+1,ny+1); Q(2:end1,2:end1) = reshape(q,nx1,ny1); clf, contourf(avg(x),avg(y),P',20,'w'), hold on contour(x,y,Q',20,'k'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k') hold off, axis equal, axis([0 lx 0 ly]) p = sort(p); caxis(p([8 end7])) title(sprintf('Re = %0.1g t = %0.2g',Re,k*dt)) drawnow end end fprintf('\n') %================================================= ====================== function B = avg(A,k) if nargin<2, k = 1; end if size(A,1)==1, A = A'; end if k<2, B = (A(2:end,:)+A(1:end1,:))/2; else, B = avg(A,k1); end if size(A,2)==1, B = B'; end function A = K1(n,h,a11) % a11: Neumann=1, Dirichlet=2, Dirichlet mid=3; A = spdiags([1 a11 0;ones(n2,1)*[1 2 1];0 a11 1],1:1,n,n)'/h^2; 
I have seen this code before, did you get an idea how to adjust it to a pipe flow problem?

I don't have alot of time on my hands, so I changed the code to inflow/outflow. You can understand what I did by using the manual written by Seibold. Check BC's and how i changed the matrices for multiplication.
Good luck! PS: It's not pipeflow, but flow between 2 parallel plates. You can check the solution by solving analytically the Navier Stokers equations for this case. You will find that the velocity is parabolic and that the pressure drop over the centerline a constant times u(0). function [Ue] = mit18086_navierstokes %MIT18086_NAVIERSTOKES % Solves the incompressible NavierStokes equations in a % rectangular domain with prescribed velocities along the % boundary. The solution method is finite differencing on % a staggered grid with implicit diffusion and a Chorin % projection method for the pressure. % Visualization is done by a colormapisoline plot for % pressure and normalized quiver and streamline plot for % the velocity field. % 07/2007 by Benjamin Seibold % http://wwwmath.mit.edu/~seibold/ % Feel free to modify for teaching and learning. % % 08/2011 by Vincent van Dijk % Changed to inflow/outflow conditions % Re = 1e0; % Reynolds number dt = 1e3; % time step tf = 1e0; % final time lx = 1; % width of box ly = 1; % height of box nx = 90; % number of xgridpoints ny = 90; % number of ygridpoints nsteps = 10; % number of steps with graphic output % nt = ceil(tf/dt); dt = tf/nt; x = linspace(0,lx,nx+1); hx = lx/nx; y = linspace(0,ly,ny+1); hy = ly/ny; [X,Y] = meshgrid(y,x); % % initial conditions U = zeros(nx1,ny); V = zeros(nx,ny1); % boundary conditions uN = x*0; vN = avg(x)*0; uS = x*0; vS = avg(x)*0; uW = avg(y)*0+1; vW = y*0; uE = avg(y)*0; vE = y*0; % Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hx^2+... [uW;zeros(nx3,ny);uE]/hy^2); Vbc = dt/Re*([vS' zeros(nx,ny3) vN']/hx^2+... [2*vW(2:end1);zeros(nx2,ny1);2*vE(2:end1)]/hy^2); fprintf('initialization') Lp = kron(speye(ny),K1(nx,hx,1,3))+kron(K1(ny,hy,1,1),s peye (nx)); perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp'; Lu = speye((nx1)*ny)+dt/Re*(kron(speye(ny),K1(nx1,hx,2,1))+... kron(K1(ny,hy,3,3),speye(nx1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru'; Lv = speye(nx*(ny1))+dt/Re*(kron(speye(ny1),K1(nx,hx,3,3))+... kron(K1(ny1,hy,2,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; %Lq = kron(speye(ny1),K1(nx1,hx,2))+kron(K1(ny1,hy,2),speye(nx1)); %perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; fprintf(', time loop\n20%%40%%60%%80%%100%%\n') for k = 1:nt % treat nonlinear terms gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ue = [2*uS'Ue(:,1) Ue 2*uN'Ue(:,end)]; Ve = [vS' V vN']; Ve = [2*vWVe(1,:);Ve;2*vEVe(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Vagamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Vagamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end1)); Ud = diff(Ue(:,2:end1))/2; Va = avg(Ve(2:end1,:)')'; Vd = diff(Ve(2:end1,:)')'/2; U2x = diff(Ua.^2gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2gamma*abs(Va).*Vd)')'/hy; U = Udt*(UVy(2:end1,:)+U2x); V = Vdt*(UVx(:,2:end1)+V2y); % implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny1); % pressure correction rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); p(perp) = Rp\(Rpt\rhs(perp)); P = reshape(p,nx,ny); U = Udiff(P)/hx; V = Vdiff(P')'/hy; uE = U(end,:);% + hx*Re*P(end1,:)/dt; Will become unstable for small dt % You can remove it because P(end,:)>0 P=P/dt; % This is the actual pressure %(notice the dt missing in the pressure correction step). % visualization if floor(25*k/nt)>floor(25*(k1)/nt), fprintf('.'), end if k==1floor(nsteps*k/nt)>floor(nsteps*(k1)/nt) % stream function %rhs = reshape(diff(U')'/hydiff(V)/hx,[],1); %q(perq) = Rq\(Rqt\rhs(perq)); %Q = zeros(nx+1,ny+1); %Q(2:end1,2:end1) = reshape(q,nx1,ny1); clf, contourf(avg(x),avg(y),P',20,'w'), hold on %contour(x,y,Q',20,'k'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k') hold off, axis equal, axis([0 lx 0 ly]) P = sort(P); caxis(P([1 end])) title(sprintf('Re = %0.1g t = %0.2g',Re,k*dt)) drawnow end end fprintf('\n') %================================================= ====================== function B = avg(A,k) if nargin<2, k = 1; end if size(A,1)==1, A = A'; end if k<2, B = (A(2:end,:)+A(1:end1,:))/2; else, B = avg(A,k1); end if size(A,2)==1, B = B'; end function A = K1(n,h,a11,a22) % a11: Neumann=1, Dirichlet=2, Dirichlet mid=3; A = spdiags([1 a11 0;ones(n2,1)*[1 2 1];0 a22 1],1:1,n,n)'/h^2; 
Hey,i am also working on SIMPLE by matlab.
I also failed to modify the cavity flow to duct flow, now I found the reason. Plz refer to my latest thread. good luck 
Thanks :)
I have noticed the issure about your second tip before. When solving for the duct flow problem, the pressure goes up to a unreasonably huge value and the solver fails to converge. So treatment has to be done to make the mass flow rate on th einlet and outlet consistent. Applying the pressure boundary instead of outflow boundary seems to require more care and attention. 
Hello, Vicent
Do you think the formulation for Ubc and Vbc has a small error? For instance, the original code is Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hx^2+[uW;zeros(nx3,ny);uE]/hy^2); but I feel the hx and hy should be switched, namely Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hy^2+[uW;zeros(nx3,ny);uE]/hx^2); Thanks. 
Fractional Step Methos
Hello Cici,
I would like to know which fractional step method are you using for your problem. I am facing the same problem as yours when I tried modifying my Lid Driven Cavity code to a duct flow. So were you able to fix the problem in your code? If yes do post the changes you had made to your code. Regards Prasanth P 
Hi Prasanth,
I was using the predictor and correct method. I have checked all the things I can think of that may result in the failure of convergence, and it turns out the boundary conditions are OK, but the linear solver may be the issue. I have written the solver in an iterative way, like for i=? j=?, not in the matrix form. I haven't verified whether this difference is responsible for the problem. If you read through the code written by an MIT professor (posted above), you will find it does work for the channel flow. Best 
Hi Cici,
My code is in FORTRAN so I have to use iterative method and I cant use the Matrix method. Regarding the MIT profs code, I assume that you implemented the out flow boundary condition and not periodic boundary condition. So did you implement that P=0 or dp/dn = 0 at the out flow along with du/dn=0 and dv/dn=0? And you did raise an issue regarding the Bcn implemented by Vincent. So is there an error? Cheers 
All times are GMT 4. The time now is 10:22. 