CFD Online Discussion Forums

CFD Online Discussion Forums (
-   Main CFD Forum (
-   -   symmetry boundary condition? (

new_at_this December 19, 2011 14:05

symmetry boundary condition?
I am solving the incompressible navier stokes equations for flow in a cylinder using matlab. My code uses a staggered grid and a fractional step method. For u I need to define a symmetric BC on the axis of symmetry so I have defined the ghost cell to equal the value of the first interior point.

However the velocity profile that it produces near the axis of symmetry is slightly less than what I expect. I think this is because I have based the BC off of the u field from the previous time step and so at every time step the ghost cell isn't exactly equal, but slightly less.

Does anyone know how to fix this issue?

houkensjtu December 19, 2011 19:53

i have exactly the same problem with symmetry boundary, hope someone can help

mettler December 20, 2011 08:50

Are you using the streamfunction-vorticity equations?

new_at_this December 20, 2011 10:51

nope, i coded in primitive variables

new_at_this December 20, 2011 13:49

i've found out that my code gets better results if I increase the number of grid point. This makes sense since if the grid points are further from the axis of symmetry, then the flux is not exactly zero between the inter point and the axis. This picture shows what I expect to get.

the profiles go from top to bottom across the axis of symmetry. When i run my code for only half the domain and apply a symmetric boundary condition, i get this...

near the axis of symmetry, the profiles start behaving irregularly which I believe is due to the boundary condition being defined from the last time step. anyone have this problem before and did you find a fix?

mettler December 20, 2011 15:35

what are you using for the symmetry BCs now?

new_at_this December 20, 2011 15:53

bottom BC: u(i-1) = u (i) or u(ghost) = u(interior), v = 0

These are the BCs that i'm using for those pics. They are the same as before just with a higher grid resolution

mettler December 20, 2011 16:25

i is the y-direction? So, your BC says that the u-velocity one step below the symmetry line equals the velocity at the symmetry line? That is not correct. Have you tried setting the symmetry-line velocity equal to Umax?

new_at_this December 20, 2011 16:40

i is in the y direction. But my ghost point below the line of symmetry is actually equal a point the same distance above the line of symmetry since it was done on a staggered grid. u doesn't exist on the symmetry line. I tried setting it equal to Umax but that produces the same irregularities.

mettler December 20, 2011 16:54

So, your symmetry line is a half-step? For example, if i = 1 then, ghost cell is at i=0 and symmetry is i=.5? Have you tried putting the symmetry line at i=0, then BCs would be Ui+1 = Ui-1, and Ui = Umax - when i=0.

new_at_this December 20, 2011 17:19

my symmetry line is at 0.5 and U is only defined at i=1,2,3,4,5.... So when the symmetry is defined at 0.5, U at 0 equals U at 1 which I think is the equivalent of what you are saying. Shifting the line of symmetry to i=0 would just change the indicies a bit, but i don't think it would affect my issue.

mettler December 20, 2011 17:26

nowhere in your code do you say that the symmetry line has a defined velocity? I think that might be your problem. I am not too sure because I have always had my symmetry line at 0, so I could define it in the code. Do you have the centerline velocity defined anywhere? It looks like your nodes close to the symmetry lines aren't getting affected by the centerline velocity like they should be..that is why they look funky..maybe???

new_at_this December 20, 2011 19:13

I don't define the center line velocity but i du/dy across the line of symmetry which from what I read should be a sufficient boundary condition. I think the issue is that I did not change the pressure boundary condition. Normally I have neumann BCs along the edges, but I don't think that is correct for the axis of symmetry and i'm not quite sure what it should be...

duri December 21, 2011 14:15


Originally Posted by new_at_this (Post 336509)
bottom BC: u(i-1) = u (i) or u(ghost) = u(interior), v = 0

These are the BCs that i'm using for those pics. They are the same as before just with a higher grid resolution

I assume you are using a Cartesian mesh. Then u(ghost)=u(interior) is correct but v(ghost)=-v(interior). All other quantities like density, pressure and enthalpy at ghost is same as interior.
For cross checking, continuity and energy flux on the symmetry face should be zero, momentum flux will be non zero due to pressure components.

mettler December 21, 2011 14:21

Also, using the CL at .5 instead of 0 has shown to cause a residual truncation error. This error is in the v-flux and v-momentum across the symmetry line. According to Roache in Computational Fluid Dynamics; "it is not preferable to specifically set the flux term equal to zero in this case."

new_at_this December 21, 2011 15:08

I think that my problem is the residual truncation error mentioned by mettler. Does Roache mention what the flux should be set to? It seems to me that it would be dependent on how far away the first interior node is from the axis of symmetry but in general du/dy seems like it should be a small number slightly greater than zero

mettler December 21, 2011 15:15

no, but he mentions that the error will be on the order of deltaY^3. I think you mentioned that if you have a larger grid that the error is not as much - that would point to the truncation error per Roach.

new_at_this December 21, 2011 16:07

Thanks, I actually went out and borrowed the book. I think there is a section that explains how to define the boundary condition for the staggered grid. I have yet to try it out though

mettler December 21, 2011 16:15

its a great book. I have a very early addition that I got from a professor when he retired..1972 edition. There is a really good section on the streamfunction-vorticity method and its related BCs. There is also a boundary condition section at the end of the book - that is where I found the information I quoted above. I used the book extensively when I did my dissertation research and still use it now.

Another good book is by Anderson (he brings in heat transfer). I got that book off of Amazon for $9.89. I honestly think that Amazon (or someone) misplaced the decimal, as the book regularly sells for closer to $100 ($98.90)

good luck

new_at_this December 22, 2011 00:57

so this is code is driving me crazy. I'm going to post it up here and if someone has some time to look at it I would really appreciate it. The code is based off the lid driven cavity written by Benjamin Seibold. It has been discussed before on the forum so I'm hoping someone is familiar with it.

I have modified the code so that both the lid and the base have velocity which produces a symmetric flow. If you copy the whole code and then run "testcode(0)" in the command line it will compute the entire domain. If you run "testcode(1)" then it will run the code for only one half of the domain using a symmetry BC on the lower boundary.


function []=testcode(sym)

close all

%% Paramters
% Input Parameters --------------------------------------------------------
AR=0.5;                                    % Aspect Ratio
Pe=500;                                      % Peclet Number
Re = 100;                                    % Reynolds number
tf = 40e-0;                                  % final time
gridnum=50;                                % Number of grid points
dt = 0.01;                                  % time step
nsteps = 10;                      % number of steps with graphic output
mode=1;                      % 1: basic , 2: extended with more plots
% Run parameters ----------------------------------------------------------

lx = 1;                                      % width of non-D domain
ly = 1;                                      % height of nond-D domain
nx=gridnum;                                  % number of x gridpoints
ny=gridnum;                                  % number of y gridpoints
hy=1/ny;                                    % Actual y grid size
hx=1/nx;                                    % Actual x grid size
% Dimensional Variables
H=1;                                      % Height of droplet
L=H/AR;                                      % Length of droplet
dy=H/ny;                                    % Width of y grid
dx=L/nx;                                    % Width of x grid
dimx(1:nx+1)=(0:nx)*dx;                      % X coordinates
dimy(1:ny+1)=(0:ny)*dy;                      % Y coordinates
% Fluid Properties
rho=1000;                                    % Density
mu=0.86e-3;                                  % Viscosity
c=4184;                                      % Specific heat
k=0.6;                                      % Thermal conductivity
ubulk=(Re*mu)/(rho*H);                      % Bulk Velocity


nt = ceil(tf/dt); dt = tf/nt;                % number of time iterations

x = linspace(0,lx,nx+1);                    % non-D x coordinates
y = linspace(0,ly,ny+1);                    % non-D y coordinates
[X,Y] = meshgrid(y,x);
% initial conditions
U = zeros(nx-1,ny); V = zeros(nx,ny-1);
% boundary conditions
uN = x*0+1; uN = [0 uN(2:end-1) 0];
if sym == 0
uS = x*0+1; uS = [0 uS(2:end-1) 0];
uW = avg(y)*0;
uE = avg(y)*0;

vN = avg(x)*0;
vS = avg(x)*0;
vW = y*0;
vE = y*0;

%% Time stepping
if sym == 0
Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+...
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+...

if sym == 1
    Lp = kron(speye(ny),AR*K1(nx,hx,1))+kron((1/AR)*K7(ny,hy,1),speye(nx));
    Lp(1,1) = 3/2*Lp(1,1);
    Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),AR^2*K1(nx-1,hx,2))+...
elseif sym == 0
    Lp = kron(speye(ny),AR*K1(nx,hx,1))+kron((1/AR)*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),AR^2*K1(nx-1,hx,2))+...

peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),AR^2*K1(nx,hx,3))+...
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
    if sym == 1
    uS = [0 U(:,1)' 0];
    Ubc = dt/Re*([0*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+...
    Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+...

    % 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*AR*(UVy(2:end-1,:)+U2x);
    V = V-dt*AR*(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);
    if sym == 0
        p(perp) = -Rp\(Rpt\rhs(perp));
        P = reshape(p,nx,ny);
    elseif sym == 1
    U = U-AR*diff(P)/hx;
    V = V-(1/AR)*diff((P)')'/hy;

    % Store U and V including ghost points
    vghost=[vW; [vS' V vN'] ;vE];
    ughost=[uS' [uW;U;uE] uN'];
    Ue = [uS' avg([uW;U;uE]')' uN'];
    Ve = [vW;avg([vS' V vN']);vE];
    % 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
        contourf(dimx,dimy,Q',20,'w-'), hold on
        Ue = [uS' avg([uW;U;uE]')' uN'];
        Ve = [vW;avg([vS' V vN']);vE];
        Len = sqrt(Ue.^2+Ve.^2+eps);
        hold off, axis equal, axis([0 L 0 H])
        %      p = sort(p); caxis(p([8 end-7]))
        title(sprintf('Re = %0.1g  t = %0.2g',Re,k*dt))




%% Profiles

axis equal, axis([0 L 0 H])
title('Droplet Streamlines (Re = 100, AR = 0.5)','Fontsize',20)

if cx1==1

if mode==2
title('U profiles Droplet Frame(Re = 100, AR = 0.5)','Fontsize',16)

title('U profiles Lab Frame (Re = 100, AR = 0.5)','Fontsize',16)

if mode==2
title('V profiles along vertical cut')

title('V profiles along horizontal cut')

%% Divergence Test
max_div = max(max(div(2:end-1,2:end-1)))

%% Vorticity Test
if mode==2
[CURLZ, CAV] = curl(X,Y,ufield,vfield);

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;

function A = K7(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 1 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;

function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = conv2(A,[1;1]/2,'valid'); else B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

All times are GMT -4. The time now is 06:56.