|
[Sponsors] |
Problem with SIMPLEC-like finite volume channel flow boundary conditions |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
June 8, 2015, 12:28 |
Problem with SIMPLEC-like finite volume channel flow boundary conditions
|
#1 |
New Member
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13 |
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; |
|
June 15, 2015, 09:26 |
|
#2 |
New Member
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13 |
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. |
|
June 15, 2015, 10:59 |
|
#3 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,839
Rep Power: 73 |
Quote:
from the plot, I see a wrong inflow ... the two boundary layers developing along the walls are not the same |
||
June 15, 2015, 11:14 |
|
#4 |
New Member
Gustavo
Join Date: Jun 2013
Posts: 26
Rep Power: 13 |
||
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |