
[Sponsors] 
August 31, 2011, 13:24 

#21 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Hi there,
I defined the boundary conditions for the outlet like this: P = 0, du/dn=0, dv/dn=0 As to the error I mentioned, I still think I am right. Regards 

September 1, 2011, 00:18 

#22 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hi,
I went through the code. I think you are right regarding the error. I guess the issue is with the Poisson solver that we are using. I am using Gauss JAcobi method. What method are you using for solving the Poisson equation? 

September 5, 2011, 10:50 

#23 
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
My apologies for not responding to this thread, for some reason I didn't get any emails
about it anymore. I would like to talk a bit about the 'boundary condition' issue cici was talking about. First of all the 'boundary condition' is only for the implicit viscosity step, it is not used anywhere else. Secondly I had already noticed that the code only works correctly for hx=hy. Though this might not be been intended, it is faster to calculate on a completely regular grid. So this potential bug does not have to cause any problems while using the code. Whats now more interesting is to see if we can derive these boundary conditions for the implicit viscosity step. So let's do that! We are calculating the matrix multiplication: (I + R(K2D))U^{n+1/2}=U^{n} Where n+1/2 denotes a time between two timesteps, K2D is the matrix approximating \frac{\partial^2}{\partial x^2} \frac{\partial^2}{\partial y^2}. I is the identity matrix (ones on the diagonal). Please note the minus sign in K2D! The 2nd order approximations are done in a centered way: U_{xx} = \frac{U_{i1,j}  2U{i,j} + U{i+1,j}}{h_x^2} and it becomes clear that boundary points are needed to calculate U_{xx} at every U location (namely 2 extra rows). Now we have to think of what the appropriate boundary conditions are. The strange thing is is that I would suggest using boundary conditions for the U field, since we are adding to this field. The equation U+Ubc where the former is a velocity and the second one is a velocity/length^2 doesn't make alot of sence to me (yet). In case of proper boundary conditions for the U field I would suggest again applying the BCs, so U=1 at inflow U = U(:,end) at outflow and U = U(inner) at the wall boundaries. It seems that Ubc is actually a boundary condition for the second derivative instead of a boundary condition for U, however I do not understand this formulation. Looking at the matrix of Ubc I notice it has the same size as U, so it does not stick boundary conditions on the outside of U but on the rim of U. Sadly I'm out of time at the moment, but I'll get back on this. I'll let you guys know when I figure it out, in the meanwhile I would love to hear your suggestions as to why it is correct or false. Regards, Vincent 

September 5, 2011, 13:12 

#24 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Hi Vincent,
I think I can answer your questions about U+Ubc. If you work through the code, you will find that Ubc is not actually for applying the boundary conditions in the normal sense, it is really for computing the RHS, or say the source term of the system of equations, which is U+Ubc. Notice that in the notes about this code, the second step is for implicit viscosity, whose formulation is (U**  U*)/dt = (Uxx** + Uyy**)/Re, simple rearrangement of terms yields U**  dt/Re*(Uxx** + Uyy**) = U*. The U* is the U in U+Ubc. How about Ubc? If we implement the above formulation to the boundary nodes, some terms need to be moved to the RHS too since some unknowns for interior nodes become known for boundary nodes. Therefore the final complete RHS for solving this ystem is U+Ubc. However, the error I described in previous thread does exist. Hope this helps. If not clear enough, I can type in the complete formulation if I've got time. Cici 

September 6, 2011, 00:10 
Regarding the code

#25 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hi Vincent,
"Secondly I had already noticed that the code only works correctly for hx=hy. Though this might not be been intended, it is faster to calculate on a completely regular grid. So this potential bug does not have to cause any problems while using the code." This bug is because of the Ubc definition. Try writing the Ubc as posted by Cici you would find that the code would start working for hx/=hy. This bug arises due to the mixup of definitions of Ubc and Vbc. Cheers PP 

September 6, 2011, 03:19 

#26  
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
Quote:
Quote:
Quote:
that expression I would be very thankfull. U** dt/Re(Uxx**+Uyy**)=U* U** dt/Re((Ui1,j**2Ui,j** +Ui+1,j**)/hx^2+(Ui,j1** 2Ui,j**+Ui,j+1**)/hy^2)=U* Now I understand that Ui1,j, Ui+1,j, Ui,j1 and Ui,j+1 are not defined at every pont. That is why we move some points the the RHS of this equation. So for instance we add dt/Re(Ui1,j)/hx^2 to the left side of our matrix U(1,:). If we correct the mistake I understand that this is the same as uW, ditto for uE. Why, however, it is in the case of hy^2 terms that a 2 appears? Thanks for the help, Regards, Vincent PS: You can improve the code by rewriting the function B: 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,k1); end if size(A,2)==1, B = B'; end You will see that this is a faster way of calculating the same thing (makes quite a difference for large matrices). PS2: wouldn't it be nice if this forum had a latex compiler? 

September 6, 2011, 05:16 
Code and Pressure Bcn

#27 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hello Vincent,
The correct expression would be: Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hy^2+... [uW;zeros(nx3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny3) vN']/hy^2+... [2*vW(2:end1);zeros(nx2,ny1);2*vE(2:end1)]/hx^2); A 2 appears not only in the hy^2 terms, it also appears in hx^2 terms. For Ubc its in hy^2 term but in Vbc it is in hx^2 terms(from the above expressions). It appears so because in x direction the x velocity is staggered from the pressure variable and in the y direction, the y velocity is staggered from pressure. I hope this clears your doubt. Now getting back to the issue with the boundary conditions for pressure. Have any of you written an iterative solver for the pressure poisson equation? Dr. Seibold's code uses a direct method to solve the the poisson equation for the pressure. I am facing issue with my iterative procedure for solving the equations. Plus, when the out flow boundary condition is partial(u)/partial(n) = 0, how can you say that P=0 or dp/dn=0 is a correct boundary condition? Cheers PP 

September 6, 2011, 05:56 

#28  
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
Quote:
In the staggered direction the terms (U(i,j+1) for instance) can be expressed as Ui,j+1 = 2uN  Ui,j We can move the 2uN to the RHS of the equation so that we do indeed get our factor 2. However the equation on the left hand side now reads (for the boundary point) Ui,j dt/Re(.../hx^2 +(Ui,j1 3Ui,j)/hy^2) Is this indeed the correct formulation? Quote:
starting from 2D NS u_t +p_x = (u_{xx} +u{yy})/Re (u^2)_x (uv)_y Assuming steady state and using partial(u)/partial(n)=0 and partial(v)/partial(y)=0 leave us with: p_x = 0 or p = constant We can choose this constant, for instance p=0 at outflow. Can you explain what issues you are facing when dealing with the iterative solver for your poisson equation? Regards, Vincent 

September 6, 2011, 06:20 
Pressure Bcn

#29 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hello Vincent,
Your formulation is correct. Regarding the boundary condition issue, so you say that you are trying to specify the boundary condition for the steady state. So this cannot be used for an evolving flow or turbulent flow. u_t=0 >steady state u_xx,(u^2)_x=0>partial(u)/partial(x)=0 But why is u_yy and (v(u)_x)=0? I have tried both P=0 and dp/dn =0. For the former case the poisson solver takes several thousand some times even lakhs depending on the domain size. For the latter case the poisson doesnt converge and doesn't diverge either. dp/dn =0 is not a correct boundary condition when the velocity boundary condition is not dirichlet. But, I dono what is the issue with my solver. Cheers PP 

September 6, 2011, 08:36 

#30  
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
Quote:
Quote:
For steady 3D flow between 2 parallel plates u=(u,0,0) we can apply NavierStokes. u_t =0 steady state u_x= 0 from NS: p_x = u_{yy}/Re (note v = 0) lhs is only function of x and rhs only of y so p_x = constant thus p(x) = c*x+p0 We can take this p0 to be equal to zero for the outflow boundary. Quote:
Could you explain a bit more on what steps you think your solver should make? Could you specify all the boundary conditions and problem you are facing? Regards, Vincent 

September 6, 2011, 11:49 

#31 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hello Vincent,
I am not interested in a steady state boundary condition. Evolving flow need not be periodic always. The boundary condition that you have given is for a very special case where v=0. Lets forget the issue about the boundary condition for the time being. Say, I use P=0 at the outflow. With my iterative Poisson solver it takes several thousand iteration or sometimes a few lakhs per time step. When I check the continuity equation at each of the pressure nodes its no where close to zero. I am not able to spot the bug in my code. So I would like to see some one solve the channel flow problem or any inflow/outflow problem with an iterative Poisson solver. The velocity is specified at all the other three boundaries. Since the velocity Bcn is a Dirichlet at these 3 boundaries, dp/dn=0 is the boundary condition for the Poisson(tats wat I assume). I am using fractional step method for solving the gov equations. Cheers PP 

September 6, 2011, 13:19 

#32 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Hi Prasanth,
What did you use for the viscosity term? If explicit schem, try change to implicit one. Cici 

September 6, 2011, 14:09 

#33 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hello Cici,
My viscous part is implicit. Even you had faced this problem while converting your Lid Driven Cavity code for a duct flow. Is your problem solved? Cheers PP 

September 6, 2011, 14:44 

#34 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Well, I can't say it is solved. Since I think both explicit and implicit viscosity will work for this problem. However, it turns out both converge now, but the solutions from explicit method are way off. So I can't say I fully understand what is going on.


September 7, 2011, 00:36 

#35 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hello Cici,
Good to hear that. Can you tell me the modifications that you made to your initial Cavity code. How are you solving your Poisson equation? Iteratively or some kind of direct Poisson solver? Cheers PP 

September 7, 2011, 13:06 

#36 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Hello,
The modifications I made are the definition of boundary conditions and the schemes for solving the viscosity. The boundary conditions are described previously, and the scheme has been changed from fully explicit to implicit, ADI in my case. As to the Poisson equation, it is solved iteratively, just make sure that the boundary conditions for pressure are represented appropriately. Best, Cici 

September 7, 2011, 15:20 

#37 
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8 
Hi Cici,
"just make sure that the boundary conditions for pressure are represented appropriately" I did not quite get your point. Can you please explain. I am using Gauss Jacobi Method to solve my Poisson iteratively. I would like to know which technique are you using. This is just to know if the limitation is due to the iterative method that I am employing. Cheers PP 

September 8, 2011, 13:19 

#38 
New Member
Join Date: Jul 2011
Posts: 25
Rep Power: 6 
Hi,
I have employed the GaussSeidel methd for solving the pressure Poisson equation. For the boundary nodes, you have to take care of the coefficients of different terms according to the boundary conditions you defined. For instance, Coeff_e = 0; Coeff_w = 1/dx^2; Coeff_n = 1/dy^2; Coeff_s = 1/dy^2; Coeff_p = 3/dx^2 2/dy^2; if the outlet boundary is defined as P = 0. Cici 

September 13, 2011, 05:22 

#39 
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
Since I have never worked with GaussSeidel iterative solver yet, I'll break the idea down and build it up from the start.
We want to solve the poisson equation: Here I added a minus sign for convinience later on. We can rewrite this as a linear system that we want to solve: Where A is our negative second difference matrix and b is our right hand side of the first equation. In order to solve this equation iteratively, we use a preconditioner matrix P. This results in the following equation: The GaussSeidel idea is to use an preconditioner P=D+L, where D is the diagonal and L is the lower triangular part of A. Let's assume we want to solve the poisson equation in 2 dimensions. This implies a matrix K2D, which can be build by K2D = kron(I,K)+kron(K,I) Where K is the matrix approximating the minus second derivative. I think it is most insightful if I change the code to use a GaussSeidel solver for the pressure, if you have any question please don't hesitate to ask. I'll post it below since my post is too long. I do have to say it's not as fast as the direct method, but you can tune down the number of iterations (i=100). Also if you want a fast method you might want to look at multigrid (for 3D applications). As to your problem, the speed of convergence depends on the eigenvalues of the matrix M=IP\A. For convergence they need to be smaller than 1, and values close to one will take a long time to dampen. The error e, can be calculated as: Good luck! Last edited by VincentD; September 14, 2011 at 03:00. 

September 13, 2011, 05:26 

#40 
New Member
Vincent
Join Date: Jul 2011
Posts: 29
Rep Power: 6 
Enjoy!
Code:
function [Ue,P] = mit18086_GS %MIT18086_NAVIERSTOKES % Solves the incompressible NavierStokes equations in a % rectangular domain with prescribed velocities along the % boundary. The solution method is finite differencing on % a staggered grid with implicit diffusion and a Chorin % projection method for the pressure. % Visualization is done by a colormapisoline plot for % pressure and normalized quiver and streamline plot for % the velocity field. % 07/2007 by Benjamin Seibold % http://wwwmath.mit.edu/~seibold/ % Feel free to modify for teaching and learning. % % 08/2011 by Vincent van Dijk % Changed to inflow/outflow conditions % Changed poisson solver to GaussSeidel % Re = 1e0; % Reynolds number dt = 1e3; % time step tf = 1e0; % final time lx = 1; % width of box ly = 1; % height of box nx = 30; % number of xgridpoints ny = 30; % number of ygridpoints nsteps = 10; % number of steps with graphic output % nt = ceil(tf/dt); dt = tf/nt; x = linspace(0,lx,nx+1); hx = lx/nx; y = linspace(0,ly,ny+1); hy = ly/ny; [X,Y] = meshgrid(y,x); % % initial conditions U = zeros(nx1,ny); V = zeros(nx,ny1); % boundary conditions uN = x*0; vN = avg(x)*0; uS = x*0; vS = avg(x)*0; uW = avg(y)*0+1; vW = y*0; uE = avg(y)*0; vE = y*0; % Ubc = dt/Re*([2*uS(2:end1)' zeros(nx1,ny2) 2*uN(2:end1)']/hy^2+... [uW;zeros(nx3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny3) vN']/hy^2+... [2*vW(2:end1);zeros(nx2,ny1);2*vE(2:end1)]/hx^2); fprintf('initialization') Lp = kron(speye(ny),K1(nx,hx,1,3))+kron(K1(ny,hy,1,1),speye(nx)); perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp'; % Changed to use GaussSeidel for the poisson equation D = diag(diag(Lp)); L = tril(Lp)D; Pr = D+L; % The preconditioner matrix M = speye(nx*ny)Pr\Lp; % The multiplication matrix Lu = speye((nx1)*ny)+dt/Re*(kron(speye(ny),K1(nx1,hx,2,1))+... kron(K1(ny,hy,3,3),speye(nx1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru'; Lv = speye(nx*(ny1))+dt/Re*(kron(speye(ny1),K1(nx,hx,3,3))+... kron(K1(ny1,hy,2,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; %Lq = kron(speye(ny1),K1(nx1,hx,2))+kron(K1(ny1,hy,2),speye(nx1)); %perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; fprintf(', time loop\n20%%40%%60%%80%%100%%\n') for k = 1:nt % 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*vWVe(1,:);Ve;2*vEVe(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Vagamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Vagamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end1)); Ud = diff(Ue(:,2:end1))/2; Va = avg(Ve(2:end1,:)')'; Vd = diff(Ve(2:end1,:)')'/2; U2x = diff(Ua.^2gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2gamma*abs(Va).*Vd)')'/hy; U = Udt*(UVy(2:end1,:)+U2x); V = Vdt*(UVx(:,2:end1)+V2y); % implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny1); % pressure correction rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); % p(perp) = Rp\(Rpt\rhs(perp)); % P = reshape(p,nx,ny); c = Pr\rhs; p = zeros(nx*ny,1); for i = 1:100; p = M*p+c; end P = reshape(p,nx,ny); U = U+diff(P)/hx; V = V+diff(P')'/hy; uE = U(end,:);% + hx*Re*P(end1,:)/dt; Will become unstable for small dt % You can remove it because P(end,:)>0 P=P/dt; % This is the actual pressure % notice the dt missing in the pressure correction step % visualization if floor(25*k/nt)>floor(25*(k1)/nt), fprintf('.'), end if k==1floor(nsteps*k/nt)>floor(nsteps*(k1)/nt) % stream function %rhs = reshape(diff(U')'/hydiff(V)/hx,[],1); %q(perq) = Rq\(Rqt\rhs(perq)); %Q = zeros(nx+1,ny+1); %Q(2:end1,2:end1) = reshape(q,nx1,ny1); clf, contourf(avg(x),avg(y),P',20,'w'), hold on %contour(x,y,Q',20,'k'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k') hold off, axis equal, axis([0 lx 0 ly]) P = sort(P); caxis(P([1 end])) title(sprintf('Re = %0.1g t = %0.2g',Re,k*dt)) drawnow end end fprintf('\n') %================================================= ====================== 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,k1); end if size(A,2)==1, B = B'; end function A = K1(n,h,a11,a22) % a11: Neumann=1, Dirichlet=2, Dirichlet mid=3; A = spdiags([1 a11 0;ones(n2,1)*[1 2 1];0 a22 1],1:1,n,n)'/h^2; 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
code for SIMPLE algorithm  2D Lid driven cavity flow problem  Collocated grid  h_amooie  OpenFOAM Programming & Development  0  June 9, 2011 11:17 
need 3D cylindrical source code for laminar flow  S. D. Ding  Main CFD Forum  0  July 23, 2002 02:01 
What is the Better Way to Do CFD?  John C. Chien  Main CFD Forum  54  April 23, 2001 08:10 
fluid flow fundas  ram  Main CFD Forum  5  June 17, 2000 21:31 
Flow visualization vs. Calculated flow patterns  Francisco Saldarriaga  Main CFD Forum  1  August 2, 1999 23:18 