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

symmetry boundary condition?

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

Reply
 
LinkBack Thread Tools Display Modes
Old   December 19, 2011, 14:05
Default symmetry boundary condition?
  #1
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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?
new_at_this is offline   Reply With Quote

Old   December 19, 2011, 19:53
Default
  #2
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
i have exactly the same problem with symmetry boundary, hope someone can help
houkensjtu is offline   Reply With Quote

Old   December 20, 2011, 08:50
Default
  #3
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
Are you using the streamfunction-vorticity equations?
mettler is offline   Reply With Quote

Old   December 20, 2011, 10:51
Default
  #4
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
nope, i coded in primitive variables
new_at_this is offline   Reply With Quote

Old   December 20, 2011, 13:49
Default
  #5
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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?
new_at_this is offline   Reply With Quote

Old   December 20, 2011, 15:35
Default
  #6
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
what are you using for the symmetry BCs now?
mettler is offline   Reply With Quote

Old   December 20, 2011, 15:53
Default
  #7
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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
new_at_this is offline   Reply With Quote

Old   December 20, 2011, 16:25
Default
  #8
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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?
mettler is offline   Reply With Quote

Old   December 20, 2011, 16:40
Default
  #9
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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.
new_at_this is offline   Reply With Quote

Old   December 20, 2011, 16:54
Default
  #10
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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.
mettler is offline   Reply With Quote

Old   December 20, 2011, 17:19
Default
  #11
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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.
new_at_this is offline   Reply With Quote

Old   December 20, 2011, 17:26
Default
  #12
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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???
mettler is offline   Reply With Quote

Old   December 20, 2011, 19:13
Default
  #13
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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...
new_at_this is offline   Reply With Quote

Old   December 21, 2011, 14:15
Post
  #14
Senior Member
 
duri
Join Date: May 2010
Posts: 130
Rep Power: 7
duri is on a distinguished road
Quote:
Originally Posted by new_at_this View Post
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.
duri is offline   Reply With Quote

Old   December 21, 2011, 14:21
Default
  #15
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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."
mettler is offline   Reply With Quote

Old   December 21, 2011, 15:08
Default
  #16
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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
new_at_this is offline   Reply With Quote

Old   December 21, 2011, 15:15
Default
  #17
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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.
mettler is offline   Reply With Quote

Old   December 21, 2011, 16:07
Default
  #18
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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
new_at_this is offline   Reply With Quote

Old   December 21, 2011, 16:15
Default
  #19
Senior Member
 
Andrew
Join Date: Mar 2009
Location: Washington, DC
Posts: 195
Rep Power: 8
mettler is on a distinguished road
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
mettler is offline   Reply With Quote

Old   December 22, 2011, 00:57
Default
  #20
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
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.

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

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Boundary Conditions Thomas P. Abraham Main CFD Forum 20 July 7, 2013 05:05
inlet velocity boundary condition murali CFX 5 August 3, 2012 08:56
Mixed/Robin boundary condition aaev OpenFOAM Bugs 2 December 15, 2011 15:03
Slip boundary condition what is inside normunds OpenFOAM Running, Solving & CFD 2 June 4, 2007 06:45
Regarding SYMMETRY boundary condition Praveen Athanki FLUENT 0 March 27, 2000 13:30


All times are GMT -4. The time now is 20:31.