[cici]

Hi All,

I have written a MATLAB code for a 2-D 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.

[VincentD]

Well done writing your own Matlab code!

If all is well (you are correctly solving Navier-Stokes 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 (free-slip) 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 free-slip 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).

[cici]

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 (non-dimensionalized), 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.

[VincentD]

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 y-direction) 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(n-1,1),1) - diag(ones(n-1,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!

[cici]

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?

[VincentD]

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

[cici]

Hi, I have found a MATLAB code online for the lid-driven 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

[VincentD]

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

[cici]

Here it is. there are 2 m-files.

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 Lid-Driven 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
%%% 288-332 - Complete B.C. for Velocity
%%% 431-435 - 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=(I-1)*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=1e-4; % 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.e-2; % 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=NI-1; % NI - 1
NJM=NJ-1; % NJ - 1
NIJ=NI*NJ; % Total number of points
%%% Indexing vector for k=(i-1)*NJ+j indexing style
LI=([1:NI]-1)*NJ; %
%%% X cell center coordinate
XC=X; XC(2:NIM)=0.5*(X(1:NIM-1)+X(2:NIM));
%%% Y cell center coordinate
YC=Y; YC(2:NJM)=0.5*(Y(1:NJM-1)+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.e-15; % 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)-----UMON-----VMON-----PMON----']);
%
%
%%% 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:NIM-1,
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(J-1); % 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 )=+CE-D; % Coefficient for E in P
AW(IJE)=-CP-D; % Coefficient for W in E
%
SU(IJ )=SU(IJ )+GDS*(FUUDS-FUCDS); % Source term for P, U
SU(IJE)=SU(IJE)-GDS*(FUUDS-FUCDS); % Source term for E, U
SV(IJ )=SV(IJ )+GDS*(FVUDS-FVCDS); % Source term for P, V
SV(IJE)=SV(IJE)-GDS*(FVUDS-FVCDS); % Source term for E, V
end;
end;
%%% FLUXES THROUGH INTERNAL NORTH CV FACES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for J=2:NJM-1,
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(I-1); % 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 )=+CN-D; % Coefficient for E
AS(IJN)=-CP-D; % Coefficient for W
%
SU(IJ )=SU(IJ )+GDS*(FUUDS-FUCDS); % Source term for P, U
SU(IJN)=SU(IJN)-GDS*(FUUDS-FUCDS); % Source term for N, U
SV(IJ )=SV(IJ )+GDS*(FVUDS-FVCDS); % Source term for P, V
SV(IJN)=SV(IJN)-GDS*(FVUDS-FVCDS); % Source term for N, V
end;
end;

%%% VOLUME INTEGRALS (SOURCE TERMS) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=2:NIM, % Loop I direction
DX=X(I)-X(I-1); % Cell length
for J=2:NJM, % Loop J direction
DY=Y(J)-Y(J-1); % 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(I-1)+P(IJ-NJ)*(1.-FX(I-1)); % Pressure at 'w'
PN=P(IJ+ 1)*FY(J )+P(IJ )*(1.-FY(J )); % Pressure at 'n'
PS=P(IJ )*FY(J-1)+P(IJ- 1)*(1.-FY(J-1)); % Pressure at 's'
DPX(IJ)=(PE-PW)/DX; % Pressure X-gradient
DPY(IJ)=(PN-PS)/DY; % Pressure Y-gradient
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=IJ-1; % 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(I-1))/(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(I-1))/(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=IJ-NJ;
U(IJB)=0.;
V(IJB)=0.;
D=VISC*(Y(J)-Y(J-1))/(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(J-1))/(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 U-VELOCITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 V-VELOCITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%

[cici]

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:NIM-1, % 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(J-1); % 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 U-velocity 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=UEL-APUE*VOLE*(DPXE-DPXEL); % 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:NJM-1, % 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(I-1); % 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 V-velocity 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=VNL-APVN*VOLN*(DPYN-DPYNL); % 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(IJ-NJ)-F1(IJ)+F2(IJ-1)-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(IJ-1)+(PP(IJ-1)-PP(IJ-2))*(1.-FY(NJM-1));
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(IJ-NJ)+(PP(IJ-NJ)-PP(IJ-NJ2))*(1.-FX(NIM-1));
end;
%%% REFERENCE PRESSURE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
IJPREF=LI(IPR)+JPR;
PP0=PP(IJPREF);
%%% CORRECT EAST MASS FLUXES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for I=2:NIM-1
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)+NJM-1,
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(I-1);
%
for J=2:NJM,
IJ=LI(I)+J;
DY=Y(J)-Y(J-1);
%
PPE=PP(IJ+NJ)*FX(I )+PP(IJ )*(1.-FX(I )); % Pressure at 'e'
PPW=PP(IJ )*FX(I-1)+PP(IJ-NJ)*(1.-FX(I-1)); % Pressure at 'w'
PPN=PP(IJ+1 )*FY(J )+PP(IJ )*(1.-FY(J )); % Pressure at 'n'
PPS=PP(IJ )*FY(J-1)+PP(IJ-1 )*(1.-FY(J-1)); % Pressure at 's'
%
U(IJ)=U(IJ)-(PPE-PPW)*DY*APR(IJ); % U-velocity corrected
V(IJ)=V(IJ)-(PPN-PPS)*DX*APR(IJ); % V-velocity 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(IJ-1)+(P(IJ-1)-P(IJ-2))*(1.-FY(NJM-1));
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(IJ-NJ)+(P(IJ-NJ)-P(IJ-NJ2))*(1.-FX(NIM-1));
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(IJ-NJ));
LS(IJ)=AS(IJ)/(1.+ALFA*UE(IJ-1));
P1=ALFA*LW(IJ)*UN(IJ-NJ);
P2=ALFA*LS(IJ)*UE(IJ-1);
LPR(IJ)=1./(AP(IJ)+P1+P2-LW(IJ)*UE(IJ-NJ)-LS(IJ)*UN(IJ-1));
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(IJ-1)-AE(IJ)*FI(IJ+NJ)...
-AW(IJ)*FI(IJ-NJ)-AP(IJ)*FI(IJ);
RESL=RESL+abs(RES(IJ));
RES(IJ)=(RES(IJ)-LS(IJ)*RES(IJ-1)-LW(IJ)*RES(IJ-NJ))*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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%

[cici]

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]);

[VincentD]

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://www-math.mit.edu/cse/codes/mi...vierstokes.pdf
The m file can be found at: http://www-math.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 Navier-Stokes 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 colormap-isoline 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://www-math.mit.edu/~seibold/
% Feel free to modify for teaching and learning.
%-----------------------------------------------------------------------
Re = 1e2; % Reynolds number
dt = 1e-2; % time step
tf = 4e-0; % final time
lx = 1; % width of box
ly = 1; % height of box
nx = 90; % number of x-gridpoints
ny = 90; % number of y-gridpoints
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(nx-1,ny); V = zeros(nx,ny-1);
% 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:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+...
[uW;zeros(nx-3,ny);uE]/hy^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+...
[2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/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((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2))+...
kron(K1(ny,hy,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3))+...
kron(K1(ny-1,hy,2),speye(nx)));
perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv';
Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1));
perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq';

fprintf(', time loop\n--20%%--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*vW-Ve(1,:);Ve;2*vE-Ve(end,:)];
Ua = avg(Ue')'; Ud = diff(Ue')'/2;
Va = avg(Ve); Vd = diff(Ve)/2;
UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx;
UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy;
Ua = avg(Ue(:,2:end-1)); Ud = diff(Ue(:,2:end-1))/2;
Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2;
U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx;
V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy;
U = U-dt*(UVy(2:end-1,:)+U2x);
V = V-dt*(UVx(:,2:end-1)+V2y);

% implicit viscosity
rhs = reshape(U+Ubc,[],1);
u(peru) = Ru\(Rut\rhs(peru));
U = reshape(u,nx-1,ny);
rhs = reshape(V+Vbc,[],1);
v(perv) = Rv\(Rvt\rhs(perv));
V = reshape(v,nx,ny-1);

% 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 = U-diff(P)/hx;
V = V-diff(P')'/hy;

% visualization
if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end
if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
% stream function
rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1);
q(perq) = Rq\(Rqt\rhs(perq));
Q = zeros(nx+1,ny+1);
Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1);
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 end-7]))
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:end-1,:))/2; else, B = avg(A,k-1); 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(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;

[cici]

I have seen this code before, did you get an idea how to adjust it to a pipe flow problem?

[VincentD]

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 Navier-Stokes 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 colormap-isoline plot for
% pressure and normalized quiver and streamline plot for
% the velocity field.

% 07/2007 by Benjamin Seibold
% http://www-math.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 = 1e-3; % time step
tf = 1e-0; % final time
lx = 1; % width of box
ly = 1; % height of box
nx = 90; % number of x-gridpoints
ny = 90; % number of y-gridpoints
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(nx-1,ny); V = zeros(nx,ny-1);
% 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:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+...
[uW;zeros(nx-3,ny);uE]/hy^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+...
[2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/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((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2,1))+...
kron(K1(ny,hy,3,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3,3))+...
kron(K1(ny-1,hy,2,2),speye(nx)));
perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv';
%Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1));
%perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq';

fprintf(', time loop\n--20%%--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*vW-Ve(1,:);Ve;2*vE-Ve(end,:)];
Ua = avg(Ue')'; Ud = diff(Ue')'/2;
Va = avg(Ve); Vd = diff(Ve)/2;
UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx;
UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy;
Ua = avg(Ue(:,2:end-1)); Ud = diff(Ue(:,2:end-1))/2;
Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2;
U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx;
V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy;
U = U-dt*(UVy(2:end-1,:)+U2x);
V = V-dt*(UVx(:,2:end-1)+V2y);

% implicit viscosity
rhs = reshape(U+Ubc,[],1);
u(peru) = Ru\(Rut\rhs(peru));
U = reshape(u,nx-1,ny);
rhs = reshape(V+Vbc,[],1);
v(perv) = Rv\(Rvt\rhs(perv));
V = reshape(v,nx,ny-1);

% 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 = U-diff(P)/hx;
V = V-diff(P')'/hy;

uE = U(end,:);% + hx*Re*P(end-1,:)/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*(k-1)/nt), fprintf('.'), end
if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
% stream function
%rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1);
%q(perq) = Rq\(Rqt\rhs(perq));
%Q = zeros(nx+1,ny+1);
%Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1);
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:end-1,:))/2; else, B = avg(A,k-1); 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(n-2,1)*[-1 2 -1];0 a22 -1],-1:1,n,n)'/h^2;

[houkensjtu]

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

[cici]

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.

[cici]

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:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+[uW;zeros(nx-3,ny);uE]/hy^2); but I feel the hx and hy should be switched, namely Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+[uW;zeros(nx-3,ny);uE]/hx^2);

Thanks.

[prasanthnitt]

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

[cici]

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

[prasanthnitt]

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