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 |
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? |
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 |
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 |
Regarding the code
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 |
Quote:
Quote:
Quote:
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? |
Code and Pressure Bcn
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 |
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,j-1 -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 |
Pressure Bcn
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 |
Quote:
Quote:
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. 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 |
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 |
Hi Prasanth,
What did you use for the viscosity term? If explicit schem, try change to implicit one. Cici |
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 |
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.
|
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 |
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 |
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 |
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 |
Enjoy!
Code:
|
All times are GMT -4. The time now is 16:59. |