[new_at_this]

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]

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

[mettler]

Are you using the streamfunction-vorticity equations?

[new_at_this]

nope, i coded in primitive variables

[new_at_this]

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.

http://s1118.photobucket.com/albums/...ent=u_prof.jpg

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...

http://s1118.photobucket.com/albums/...nt=symprof.jpg

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]

what are you using for the symmetry BCs now?

[new_at_this]

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]

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]

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]

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]

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]

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]

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]

Originally Posted by new_at_this

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]

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]

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]

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]

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]

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]

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
clc

%% 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];
end
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+...
    [uW;zeros(nx-3,ny);uE]/hx^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+...
    [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2);
end

fprintf('initialization')
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))+...
    kron(K7(ny,hy,3),speye(nx-1)));
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))+...
    kron(K1(ny,hy,3),speye(nx-1)));
end

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))+...
    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
    
    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+...
        [uW;zeros(nx-3,ny);uE]/hx^2);
    Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+...
        [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2);
    end

    % 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
        p2=Lp\rhs;
        P2=reshape(p2,nx,ny);
        P=-P2;
    end
    
    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);
        quiver(dimx,dimy,(Ue./Len)',(Ve./Len)',.4,'k-')
        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))
        drawnow

    end
end
fprintf('\n')

ufield=Ue';
vfield=Ve';

ulab=abs(ufield-1);

%% Profiles

figure(1)
contour(dimx,dimy,Q',20,'Linewidth',1.2);
axis equal, axis([0 L 0 H])
set(gca,'FontSize',14)
xlabel('X','Fontsize',16)
ylabel('Y','Fontsize',16)
title('Droplet Streamlines (Re = 100, AR = 0.5)','Fontsize',20)
colorbar



mid=ceil((nx+1)/2);
cx1=ceil(1*(nx+1)/100);
if cx1==1
    cx1=2;
end
cx2=ceil(5*(nx+1)/100);
cx3=ceil(1*(nx+1)/5);

if mode==2
figure(2)
plot(ufield(:,cx1),y,'k-',ufield(:,cx2),y,':',ufield(:,cx3),y,'m-.',ufield(:,mid),y,'--');
leg=legend('0.02H','0.1H','0.4L','H','other','location','East');
xlabel('U^*','Fontsize',16)
ylabel('Y','Fontsize',16)
title('U profiles Droplet Frame(Re = 100, AR = 0.5)','Fontsize',16)
set(leg,'FontSize',14);
end

figure(3)
uprof=plot(ulab(:,cx1),y,'k-',ulab(:,cx2),y,':',ulab(:,cx3),y,'m-.',ulab(:,mid),y,'--');
set(uprof,'LineWidth',1.5)
set(gca,'FontSize',14)
leg=legend('0.02H','0.1H','0.4L','H','other','location','West');
xlabel('U^*','Fontsize',16)
ylabel('Y','Fontsize',16)
title('U profiles Lab Frame (Re = 100, AR = 0.5)','Fontsize',16)
set(leg,'FontSize',14);

if mode==2
figure(4)
plot(vfield(:,cx1),x,vfield(:,cx2),x,vfield(:,cx3),x,vfield(:,mid),x)%,vfield(:,end-quart),y)
legend('0.02H','0.1H','0.4L','H','other','location','NorthWest')
title('V profiles along vertical cut')

figure(5)
plot(x,vfield(cx1,:),x,vfield(cx2,:),x,vfield(cx3,:),x,vfield(mid,:))%,vfield(:,end-quart),y)
legend('0.02H','0.1H','0.4L','H','other','location','NorthWest')
title('V profiles along horizontal cut')
end


%% Divergence Test
div=divergence(ufield,vfield);
max_div = max(max(div(2:end-1,2:end-1)))

%% Vorticity Test
if mode==2
[CURLZ, CAV] = curl(X,Y,ufield,vfield);
figure(6)
contourf(x,y,CURLZ,50)
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;

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