CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   symmetry boundary condition? (http://www.cfd-online.com/Forums/main/95473-symmetry-boundary-condition.html)

 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.

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

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

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

All times are GMT -4. The time now is 01:13.