[cici]

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

[prasanthnitt]

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?

[VincentD]

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_{i-1,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

[cici]

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

[prasanthnitt]

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

[VincentD]

Originally Posted by prasanthnitt

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

Well, that is a really clear confirmation of the error made defining Ubc and Vbc. Thanks for figuring that one out.

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.

This part was completely clear, I already formulated this in my own post ("the matrix multiplication is (I +R(K2D))U^{**}=U^{*}).

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

I can't say I understand it yet, if you could derive why Ubc has to have
that expression I would be very thankfull.

U** -dt/Re(Uxx**+Uyy**)=U*
U** -dt/Re((Ui-1,j**-2Ui,j** +Ui+1,j**)/hx^2+(Ui,j-1** -2Ui,j**+Ui,j+1**)/hy^2)=U*

Now I understand that Ui-1,j, Ui+1,j, Ui,j-1 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(Ui-1,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,k-1); 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?

[prasanthnitt]

Hello Vincent,

The correct expression would be:

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

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

[VincentD]

Hello Vincent,

The correct expression would be:

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

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.

Alright I understand that we defined U and V velocity fields on a staggered grid. In the nonstaggered direction we do not have a factor, as we would expect.
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,j-1 -3Ui,j)/hy^2)

Is this indeed the correct formulation?

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

When you have partial(u)/partial(n) = 0 this also implies partial(v)/partial(y)=0 (continuity).
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

[prasanthnitt]

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

[VincentD]

Originally Posted by prasanthnitt

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.

Indeed for a turbulent flow we do not expect p=0 at the outflow (this would imply a steady field). I guess an evolving flow means periodic boundary conditions, so p is not equal to zero in that case.

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 guess you mean to ask why u_yy and vu_y are equal to zero. In my previous post I took for granted that partial(u)/partial(n) =0.

For steady 3D flow between 2 parallel plates u=(u,0,0) we can apply Navier-Stokes.

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.

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

If you prescribe \partial p / \partial n at all the boundaries then you obtain a singular matrix with a determinat of zero.
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

[prasanthnitt]

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

[cici]

Hi Prasanth,

What did you use for the viscosity term? If explicit schem, try change to implicit one.

Cici

[prasanthnitt]

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

[cici]

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.

[prasanthnitt]

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

[cici]

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

[prasanthnitt]

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

[cici]

Hi,

I have employed the Gauss-Seidel 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

[VincentD]

Since I have never worked with Gauss-Seidel 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 Gauss-Seidel 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 Gauss-Seidel 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=I-P\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!

[VincentD]

Enjoy!

 
function [Ue,P] = mit18086_GS
%MIT18086_NAVIERSTOKES
% Solves the incompressible Navier-Stokes 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 colormap-isoline plot for
% pressure and normalized quiver and streamline plot for
% the velocity field.
 
% 07/2007 by Benjamin Seibold
% http://www-math.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 Gauss-Seidel
%-----------------------------------------------------------------------
Re = 1e0; % Reynolds number
dt = 1e-3; % time step
tf = 1e-0; % final time
lx = 1; % width of box
ly = 1; % height of box
nx = 30; % number of x-gridpoints
ny = 30; % number of y-gridpoints
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(nx-1,ny); V = zeros(nx,ny-1);
% 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: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);
 
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 Gauss-Seidel 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((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2,1))+...
kron(K1(ny,hy,3,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3,3))+...
kron(K1(ny-1,hy,2,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
% 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*(UVy(2:end-1,:)+U2x);
V = V-dt*(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);
% 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(end-1,:)/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*(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
%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,k-1); 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(n-2,1)*[-1 2 -1];0 a22 -1],-1:1,n,n)'/h^2;