[ghobold]

I'm trying to implement a CFD code (Matlab) to solve a channel flow with inlet and outlet boundary conditions, but I am having problems with the outflow boundary condition (I believe). I was able to make it work for the lid-driven cavity, but I can't make it work for channel flow.

The scheme I am using is PRIME (PRessure Implicit Momentum Explicit), which is basically a modification of the SIMPLEC algorithm, which solves the momentum equations explicitly and the Poisson equation for pressure (not pressure correction) implicitly in a staggered grid. It works fine for the lid-driven cavity, so I believe the problem is really in the boundary conditions.

Here are the boundary conditions I am using:

Inlet:
u = 1, v = 0

Wall:
u = 0, v= 0

Outflow:
u(end) = u(end-1), v(end) = v(end-1), and set p(end, end) = 0;

Although the code converges quickly, and the fully developed profile is obtained at outlet, the entrance length is severely underestimated and the flow behaves as if fully developed a few volumes from the inlet.


Here's how the solution looks like, the contour is u-velocity along the channel. The entrance length is supposed to be ~0.6 and it's much smaller. Also note the anomality near the inlet

Any idea on how to make my code work? Is this a common error? Am I doing something wrong with my boundary conditions?

Here's the code I am running:

clear all
close all

H = 0.05; %channel height
L = 1; %channel length
nx = 50; %number of x volumes
ny = 50; %number of y volumes

dy = H/ny;
dx = L/nx;

x_aux = 0:dx:L;
y_aux = 0:dy:H;

[xb, yb] = meshgrid(x_aux, y_aux); yb = flipud(yb); %create mesh of border position

xp = 0.5*(xb(:, 1:end-1) + xb(:, 2:end)); xp = xp(1:end-1, :); %central position of the non-staggered volume
yp = 0.5*(yb(1:end-1, :) + yb(2:end, :)); yp = yp(:, 1:end-1); %central position of the non-staggered volume

xu = xb(1:end-1, 2:end); %central position of the x-staggered volume
yu = yp; %central position of the x-staggered volume

xv = xp; %central position of the y-staggered volume
yv = yb(2:end, 1:end-1); %central position of the y-staggered volume

% Initial guess
u = 0.01*ones(ny, nx);
v = 0.01*ones(ny, nx);
pressure = zeros(ny, nx);

% Flow parameters
Re = 100;
rho = 1;
mu = 1/Re;
U = 1;

for k = 1:1000
    %% solve for x
    ux_e = [0.5*(u(:, 2:end) + u(:, 1:end-1)), 0.5*(u(:, end) + u(:, end-1))];
    ux_w = 0.5*(u(:, 1:end) + [U*ones(ny, 1), u(:, 1:end-1)]);
    vx_n = [0.5*( v(:, 1:end-1) + v(:, 2:end) ), zeros(ny, 1)];
    vx_s = [[0.5*( v(2:end, 1:end-1) + v(2:end, 2:end) ), zeros(ny-1, 1)]; zeros(1, nx)];
    
    % Calculate advective coefficient
    M_Xe = rho .* ux_e .* dy;
    M_Xw = rho .* ux_w .* dy;
    M_Xn = rho .* vx_n .* dx;
    M_Xs = rho .* vx_s .* dx;
    
    % Calculate diffusive coefficient
    D_Xe = mu .* dy ./ dx;
    D_Xw = mu .* dy ./ dx;
    D_Xn = mu .* dx ./ dy;
    D_Xs = mu .* dx ./ dy;
    
    % Calculate volume Peclet number
    Pe_e = dx .* ux_e .* rho ./ mu;
    Pe_w = dx .* ux_w .* rho ./ mu;
    Pe_n = dy .* vx_n .* rho ./ mu;
    Pe_s = dy .* vx_s .* rho ./ mu;
    
    % Interpolation scheme (sort of a power law)
    alpha_e = sign(Pe_e).*Pe_e.^2./(10+2*Pe_e.^2);
    alpha_w = sign(Pe_w).*Pe_w.^2./(10+2*Pe_w.^2);
    alpha_n = sign(Pe_n).*Pe_n.^2./(10+2*Pe_n.^2);
    alpha_s = sign(Pe_s).*Pe_s.^2./(10+2*Pe_s.^2);
    
    beta_e = (1+0.005*Pe_e.^2)./(1+0.05*Pe_e.^2);
    beta_w = (1+0.005*Pe_w.^2)./(1+0.05*Pe_w.^2);
    beta_n = (1+0.005*Pe_n.^2)./(1+0.05*Pe_n.^2);
    beta_s = (1+0.005*Pe_s.^2)./(1+0.05*Pe_s.^2);
    
    % Calculate coefficients
    A_Xe = - (1/2  - alpha_e) .* M_Xe + D_Xe .* beta_e;
    A_Xw =  (1/2  + alpha_w) .* M_Xw + D_Xw .* beta_w;
    A_Xn = - (1/2  - alpha_n) .* M_Xn + D_Xn .* beta_n;
    A_Xs =  (1/2  + alpha_s) .* M_Xs + D_Xs .* beta_s;
    A_Xp = A_Xe + A_Xw + A_Xn + A_Xs;
    B_X = zeros(ny, nx);
    
    % Apply boundary conditions
    % left boundary u=U
    A_Xp(:, 1) = A_Xe(:, 1) + A_Xw(:, 1) + A_Xn(:, 1) + A_Xs(:, 1);
    B_X(:, 1) = B_X(:, 1) +  A_Xw(:, 1)*U;
    A_Xw(:, 1) = 0;
    
    %lower boundary u=0
    A_Xs(end, :) =  (1/2  + alpha_n(end, :)) .* M_Xs(end, :) + 2*D_Xs(end, :) .* beta_n(end, :);
    A_Xp(end, :) = A_Xe(end, :) + A_Xw(end, :) + A_Xn(end, :) + A_Xs(end, :);
    B_X(end, :) = B_X(end, :) +  A_Xs(end, :)*0;
    A_Xs(end, :) = 0;
    
    %upper boundary u=0
    A_Xn(1, :) =  - (1/2  - alpha_n(1, :)) .* M_Xn(1, :) + 2*D_Xn(1, :) .* beta_n(1, :);
    A_Xp(1, :) = A_Xe(1, :) + A_Xw(1, :) + A_Xn(1, :) + A_Xs(1, :);
    B_X(1, :) = B_X(1, :) +  A_Xn(1, :)*0;
    A_Xn(1, :) = 0;
    
    %right boundary u(end)=u(end-1)
    A_Xp(:, end) = 1;
    A_Xe(:, end) = 0;
    A_Xw(:, end) = 1;
    A_Xn(:, end) = 0;
    A_Xs(:, end) = 0;
    B_X(:, end) = 0;
    
    %north-west corner u_west=U, u_north=0
    A_Xn(1, 1) =  - (1/2  - alpha_n(1, 1)) .* M_Xn(1, 1) + 2*D_Xn(1, 1) .* beta_n(1, 1);
    A_Xp(1, 1) = A_Xe(1, 1) + A_Xw(1, 1) + A_Xn(1, 1) + A_Xs(1, 1);
    B_X(1, 1) = B_X(1, 1) +  A_Xw(1, 1)*U + A_Xn(1, 1)*0;
    A_Xn(1, 1) = 0;
    A_Xw(1, 1) = 0;
    
    %south-west corner u_west=U, u_south=0
    A_Xn(end, 1) =  - (1/2  - alpha_n(end, 1)) .* M_Xn(end, 1) + 2*D_Xn(end, 1) .* beta_n(end, 1);
    A_Xp(end, 1) = A_Xe(end, 1) + A_Xw(end, 1) + A_Xn(end, 1) + A_Xs(end, 1);
    B_X(end, 1) = B_X(end, 1) +  A_Xw(end, 1)*U + A_Xs(end, 1)*0;
    A_Xs(end, 1) = 0;
    A_Xw(end, 1) = 0;
    
    
    % Explicitly solve the x-momentum equation
    u_Xe = [u(:, 2:end), u(:, end)];
    u_Xw = [U*ones(ny, 1), u(:, 1:end-1)];
    u_Xn = [ones(1, nx)*0; u(1:end-1, :)];
    u_Xs = [u(2:end, :); zeros(1, nx)*0];
    
    uhat = (A_Xe .* u_Xe + A_Xw .* u_Xw + A_Xn .* u_Xn + A_Xs .* u_Xs + B_X)./A_Xp;
    
    %% solve for y
    
    % Define face velocities for the y-staggered grid
    uy_e = [ones(1, nx)*0; 0.5*(u(2:end, :) + u(1:end-1, :))];
    uy_w = [ones(1, nx)*0; U*ones(ny-1, 1), 0.5*(u(2:end, 1:end-1) + u(1:end-1, 1:end-1))];
    vy_n = [zeros(1, nx); 0.5*(v(1:end-1, :) + v(2:end, :))];
    vy_s = 0.5*(v(1:end, :) + [v(2:end, :);zeros(1, nx)]);
    
    % Calculate advective component
    M_Ye = rho .* uy_e .* dy;
    M_Yw = rho .* uy_w.* dy;
    M_Yn = rho .* vy_n .* dx;
    M_Ys = rho .* vy_s .* dx;
    
    % Calculate diffusive component
    D_Ye = mu .* dy ./ dx;
    D_Yw = mu .* dy ./ dx;
    D_Yn = mu .* dx ./ dy;
    D_Ys = mu .* dx ./ dy;
    
    % Calculate volume Peclet number
    Pe_e = dx .* uy_e .* rho ./ mu;
    Pe_w = dx .* uy_w .* rho ./ mu;
    Pe_n = dy .* vy_n .* rho ./ mu;
    Pe_s = dy .* vy_s .* rho ./ mu;
    
    % Interpolation scheme (sort of a power law)
    alpha_e = sign(Pe_e).*Pe_e.^2./(10+2*Pe_e.^2);
    alpha_w = sign(Pe_w).*Pe_w.^2./(10+2*Pe_w.^2);
    alpha_n = sign(Pe_n).*Pe_n.^2./(10+2*Pe_n.^2);
    alpha_s = sign(Pe_s).*Pe_s.^2./(10+2*Pe_s.^2);
    
    beta_e = (1+0.005*Pe_e.^2)./(1+0.05*Pe_e.^2);
    beta_w = (1+0.005*Pe_w.^2)./(1+0.05*Pe_w.^2);
    beta_n = (1+0.005*Pe_n.^2)./(1+0.05*Pe_n.^2);
    beta_s = (1+0.005*Pe_s.^2)./(1+0.05*Pe_s.^2);
    
    % Calculate coefficients
    A_Ye = - (1/2  - alpha_e) .* M_Ye + D_Ye .* beta_e;
    A_Yw =  (1/2  + alpha_w) .* M_Yw + D_Yw .* beta_w;
    A_Yn = - (1/2  - alpha_n) .* M_Yn + D_Yn .* beta_n;
    A_Ys =  (1/2  + alpha_s) .* M_Ys + D_Ys .* beta_s;
    A_Yp = A_Ye + A_Yw + A_Yn + A_Ys;
    B_Y = zeros(ny, nx);
    
    % Apply boundary conditions
    % left boundary v=0
    A_Yw(:, 1) =  (1/2  + alpha_n(:, 1)) .* M_Yw(:, 1) + 2*D_Yw(:, 1) .* beta_n(:, 1);
    A_Yp(:, 1) = A_Ye(:, 1) + A_Yw(:, 1) + A_Yn(:, 1) + A_Ys(:, 1);
    B_Y(:, 1) = B_Y(:, 1) +  A_Yw(:, 1)*0;
    A_Yw(:, 1) = 0;
    
    % right boundary v(end) = v(end-1)
    A_Yp(:, end) = 1;
    A_Ye(:, end) = 0;
    A_Yw(:, end) = 1;
    A_Yn(:, end) = 0;
    A_Ys(:, end) = 0;
    B_Y(:, end) = 0;
    
    %lower boundary v=0
    A_Ys(end, :) = 0;
    
    %upper boundary v=0
    A_Yp(1, :) = 1;
    A_Ye(1, :) = 0;
    A_Yw(1, :) = 0;
    A_Yn(1, :) = 0;
    A_Ys(1, :) = 0;
    B_Y(1, :) = 0;
    
    %south-west corner v_west = v_south = 0
    A_Yw(end, 1) =  - (1/2  - alpha_n(end, 1)) .* M_Yw(end, 1) + 2*D_Yw(end, 1) .* beta_n(end, 1);
    A_Yp(end, 1) = A_Ye(end, 1) + A_Yw(end, 1) + A_Yn(end, 1) + A_Ys(end, 1);
    B_Y(end, 1) = B_Y(end, 1) +  A_Ys(end, 1)*0;
    A_Ys(end, 1) = 0;
    A_Yw(end, 1) = 0;
    
    %south-east corner v_west = v_north = 0
    A_Ye(end, end) =  - (1/2  - alpha_e(end, end)) .* M_Ye(end, end) + 2*D_Ye(end, end) .* beta_e(end, end);
    A_Yp(end, end) = A_Ye(end, end) + A_Yw(end, end) + A_Yn(end, end) + A_Ys(end, end);
    B_Y(end, end) = B_Y(end, end) +  A_Ye(end, end)*0 + A_Ys(end, end)*0;
    A_Ys(end, end) = 0;
    A_Ye(end, end) = 0;
    
    v_Ye = [v(:, 2:end), v(:, end)];
    v_Yw = [zeros(ny, 1), v(:, 1:end-1)];
    v_Yn = [zeros(1, nx); v(1:end-1, :)];
    v_Ys = [v(2:end, :); zeros(1, nx)];
    
    vhat = (A_Ye .* v_Ye + A_Yw .* v_Yw + A_Yn .* v_Yn + A_Ys .* v_Ys + B_Y)./A_Yp;
    
    %% solve for continuity
    
    % Calculate d coefficients (as in SIMPLE)
    d_x = dy./A_Xp;
    d_y = dx./A_Yp;
    d_xe = d_x;
    d_xw = [zeros(ny, 1), d_x(:, 1:end-1)];
    d_yn = d_y;
    d_ys = [d_y(2:end, :); zeros(1, nx)];
    
    % Calculate the coefficients of the continuity equation
    A_Ce = -dy .* d_xe;
    A_Cw = -dy .* d_xw;
    A_Cn = -dx .* d_yn;
    A_Cs = -dx .* d_ys;
    A_Ce(:, end) = 0;
    A_Cw(:, 1) = 0;
    A_Cn(1, :) = 0;
    A_Cs(end, :) = 0;
    A_Cp = -(A_Ce + A_Cw + A_Cn + A_Cs);
    
    % Calculate the face velocities of the non-staggered mesh
    ue = uhat;
    uw = [U*ones(ny, 1), uhat(:, 1:end-1)];
    vn = vhat;
    vs = [vhat(2:end, :); zeros(1, nx)];
    
    % Calculate the velocity divergent
    B_C = -(dy*(ue - uw) + dx*(vn-vs));
    
    % Set reference pressure
    A_Ce(end, end) = 0;
    A_Cw(end, end) = 0;
    A_Cn(end, end) = 0;
    A_Cs(end, end) = 0;
    A_Cp(end, end) = 1;
    B_C(end, end) = 0;
    
    
    
    % The following piece of code concatenates the coefficients into a
    % sparse matrix in order to follow the continuity equation
    [n, m] = size(A_Cp);
    
    % Upper-left cell
    i=1;
    j=1;
    ij=i+n.*(j-1);
    row=[ij; ij; ij];
    col=[ij; ij+1; ij+n];
    a_S=reshape(A_Cs(ij), 1, 1);
    a_E=reshape(A_Ce(ij), 1, 1);
    a_C=reshape(A_Cp(ij), 1, 1);
    val=[a_C; a_S; a_E];
    
    % Bottom-left cell
    i=n;
    j=1;
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij];
    col=[col; ij; ij-1; ij+n];
    a_N=reshape(A_Cn(ij), 1, 1);
    a_E=reshape(A_Ce(ij), 1, 1);
    a_C=reshape(A_Cp(ij), 1, 1);
    val=[val; a_C; a_N; a_E];
    
    % Upper-right cell
    i=1;
    j=m;
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij];
    col=[col; ij; ij+1; ij-n];
    a_S=reshape(A_Cs(ij), 1, 1);
    a_W=reshape(A_Cw(ij), 1, 1);
    a_C=reshape(A_Cp(ij), 1, 1);
    val=[val; a_C; a_S; a_W];
    
    % Bottom-right cell
    i=n;
    j=m;
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij];
    col=[col; ij; ij-1; ij-n];
    a_N=reshape(A_Cn(ij), 1, 1);
    a_W=reshape(A_Cw(ij), 1, 1);
    a_C=reshape(A_Cp(ij), 1, 1);
    val=[val; a_C; a_N; a_W];
    
    % Upper boundary
    i=1;
    j=(2:m-1)';
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij; ij];
    col=[col; ij; ij+1;  ij-n; ij+n];
    a_S=reshape(A_Cs(ij), (m-2), 1);
    a_E=reshape(A_Ce(ij), (m-2), 1);
    a_W=reshape(A_Cw(ij), (m-2), 1);
    a_C=reshape(A_Cp(ij), (m-2), 1);
    val=[val; a_C; a_S; a_W; a_E];
    
    % Lower boundary
    i=n;
    j=(2:m-1)';
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij; ij];
    col=[col; ij; ij-1;  ij-n; ij+n];
    a_N=reshape(A_Cn(ij), (m-2), 1);
    a_E=reshape(A_Ce(ij), (m-2), 1);
    a_W=reshape(A_Cw(ij), (m-2), 1);
    a_C=reshape(A_Cp(ij), (m-2), 1);
    val=[val; a_C; a_N; a_W; a_E];
    
    % Left boundary
    i=(2:n-1)';
    j=1;
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij; ij];
    col=[col; ij; ij-1;  ij+1; ij+n];
    a_N=reshape(A_Cn(ij), (n-2), 1);
    a_S=reshape(A_Cs(ij), (n-2), 1);
    a_E=reshape(A_Ce(ij), (n-2), 1);
    a_C=reshape(A_Cp(ij), (n-2), 1);
    val=[val; a_C; a_N; a_S; a_E];
    
    % Right boundary
    i=(2:n-1)';
    j=m;
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij; ij];
    col=[col; ij; ij-1;  ij+1; ij-n];
    a_N=reshape(A_Cn(ij), (n-2), 1);
    a_S=reshape(A_Cs(ij), (n-2), 1);
    a_W=reshape(A_Cw(ij), (n-2), 1);
    a_C=reshape(A_Cp(ij), (n-2), 1);
    val=[val; a_C; a_N; a_S; a_W];
    
    % Central cells
    i=repmat((2:n-1)',m-2,1);
    j=reshape(repmat((2:m-1),n-2,1),(n-2)*(m-2),1);
    ij=i+n.*(j-1);
    row=[row; ij; ij; ij; ij; ij];
    col=[col; ij; ij-1; ij+1; ij-n; ij+n];
    a_N=reshape(A_Cn(ij), (n-2)*(m-2), 1);
    a_S=reshape(A_Cs(ij), (n-2)*(m-2), 1);
    a_E=reshape(A_Ce(ij), (n-2)*(m-2), 1);
    a_W=reshape(A_Cw(ij), (n-2)*(m-2), 1);
    a_C=reshape(A_Cp(ij), (n-2)*(m-2), 1);
    val=[val; a_C; a_N; a_S; a_W; a_E];
    
    % Define A as a sparse matrix (saves a lot of memory)
    A = sparse(row, col, val);
    
    % Reshape the right hand side matrix so that it becomes a vector
    b = reshape(B_C, n*m, 1);
    
    % Store last known pressure
    pressure_ = pressure;
    
    % Solve pressure equation
    pressure = A\b;
    pressure = reshape(pressure, n, m);
    
    % Set border pressures for velocity correction
    pE = [pressure(:, 2:end), pressure(:, end)];
    pN = [pressure(1, :); pressure(1:end-1, :)];
    
    % Store last known velocities
    u_ = u;
    v_ = v;
    
    % Velocity correction
    u = uhat + d_x .* (pressure - pE);
    v = vhat + d_y .* (pressure - pN);
    
    % Ensure global mass conservation
    m_in = H*U;
    m_out = sum(u(:, end).*dy);
    m_r = m_in - m_out;
    beta = m_r./m_out;
    u(:, end) = u(:, end)*(1+beta);
    
    % Calculate residual
    res = sum(sum(abs(u-u_) + abs(v-v_) + sum(sum(abs(abs(pressure) - abs(pressure_))))))./nx./ny
end
%% plot
contourf(xb(1:end-1, 2:end), yp, u, 100, 'linewidth', 0); colorbar;