CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Problem with SIMPLEC-like finite volume channel flow boundary conditions

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 8, 2015, 12:28
Default Problem with SIMPLEC-like finite volume channel flow boundary conditions
  #1
New Member
 
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13
ghobold is on a distinguished road
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:

Code:
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;
ghobold is offline   Reply With Quote

Old   June 15, 2015, 09:26
Default
  #2
New Member
 
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13
ghobold is on a distinguished road
For those interested in what's wrong with the code, there was nothing. I was just comparing the results with the wrong data.

The Reynolds number specified in the code is not really the flow Reynolds number, its an order lower. For example, if I specify Re = 1000, then mu = 1e-4, rho = 1, U = 1, Dh = 0.1, and thus the flow Nusselt number is 100, not 1000.

The code solves the Navier-Stokes equation in a parallel plate with inflow and outflow boundary conditions. I am now going to implement specified pressure boundary conditions.
ghobold is offline   Reply With Quote

Old   June 15, 2015, 10:59
Default
  #3
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,839
Rep Power: 73
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by ghobold View Post
For those interested in what's wrong with the code, there was nothing. I was just comparing the results with the wrong data.

The Reynolds number specified in the code is not really the flow Reynolds number, its an order lower. For example, if I specify Re = 1000, then mu = 1e-4, rho = 1, U = 1, Dh = 0.1, and thus the flow Nusselt number is 100, not 1000.

The code solves the Navier-Stokes equation in a parallel plate with inflow and outflow boundary conditions. I am now going to implement specified pressure boundary conditions.

from the plot, I see a wrong inflow ... the two boundary layers developing along the walls are not the same
FMDenaro is offline   Reply With Quote

Old   June 15, 2015, 11:14
Default
  #4
New Member
 
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13
ghobold is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
from the plot, I see a wrong inflow ... the two boundary layers developing along the walls are not the same
The inflow problem is solved if the number of iterations increases and/or mesh is refined.
ghobold is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
how to set periodic boundary conditions Ganesh FLUENT 15 November 18, 2020 06:09
Overflow Error in Multiphase Modelling with Two Continuous Fluids ashtonJ CFX 6 August 11, 2014 14:32
Radiation interface hinca CFX 15 January 26, 2014 17:11
[blockMesh] BlockMesh FOAM warning gaottino OpenFOAM Meshing & Mesh Conversion 7 July 19, 2010 14:11
[blockMesh] Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Meshing & Mesh Conversion 10 April 2, 2007 14:00


All times are GMT -4. The time now is 00:39.