|
[Sponsors] |
pressure boundary conditions for collocated grid for Navier-Stokes |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
May 17, 2017, 14:52 |
Pressure BC
|
#21 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Hi Dr. Denaro,
Is that not what I'm doing in terms of the code above? In terms of setting u* and v* I am aware there is an issue of accuracy, and to do it correctly I am aware I have to relate u* and P. Code:
for j=3:ny-2 for i=3:nx-2 R(i,j) = (rho/dt)*((1/dx)*(ufs(i,j) - ufs(i-1,j)) + (1/dy)*(vfs(i,j) - vfs(i,j-1))); end end Iter = 1; IterCntr = 0; while (IterCntr == 0 && Iter <= IterMax) for j=3:ny-2 for i=3:nx-2 P(i,j) = (1-omega)*P(i,j) + omega/(2/dx^2 + 2/dy^2)*((P(i+1,j) + P(i-1,j))/dx^2 + (P(i,j+1) + P(i,j-1))/dy^2 - R(i,j)); end end for j=3:ny-2 R(2,j) = (rho/dt)*((1/dx)*(ufs(2,j) - u(1,j)) + (1/dy)*(vfs(2,j) - vfs(2,j-1))); P(2,j) = (1-omega)*P(2,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(3,j) + (1/dy^2)*(P(2,j+1) + P(2,j-1)) - R(2,j)); R(nx-1,j) = (rho/dt)*((1/dx)*(u(nx,j) - ufs(nx-2,j)) + (1/dy)*(vfs(nx-1,j) - vfs(nx-1,j-1))); P(nx-1,j) = (1-omega)*P(nx-1,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(nx-2,j) + (1/dy^2)*(P(nx-1,j+1) + P(nx-1,j-1)) - R(nx-1,j)); end for i=3:nx-2 R(i,2) = (rho/dt)*((1/dx)*(ufs(i,2) - ufs(i-1,2)) + (1/dy)*(vfs(i,2) - v(i,1))); P(i,2) = (1-omega)*P(i,2) + omega/(2/dx^2 + 1/dy^2)*((1/dx^2)*(P(i+1,2) + P(i-1,2)) + (1/dy^2)*P(i,3) - R(i,2)); R(i,ny-1) = (rho/dt)*((1/dx)*(ufs(i,ny-1) - ufs(i-1,ny-1)) + (1/dy)*(v(i,ny) - vfs(i,ny-2))); P(i,ny-1) = (1-omega)*P(i,ny-1) + omega/(2/dx^2 + 1/dy^2)*((1/dx^2)*(P(i+1,ny-1) + P(i-1,ny-1)) + (1/dy^2)*P(i,ny-2) - R(i,ny-1)); end R(2,ny-1) = (rho/dt)*((1/dx)*(ufs(2,ny-1) - u(1,ny-1)) + (1/dy)*(v(1,ny) - vfs(2,ny-1))); P(2,ny-1) = (1-omega)*P(2,ny-1) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(3,ny-1) + (1/dy^2)*P(2,ny-2) - R(2,ny-1)); R(nx-1,2) = (rho/dt)*((1/dx)*(u(nx,2) - ufs(nx-2,2)) + (1/dy)*(vfs(nx-1,2) - v(nx-1,1))); P(nx-1,2) = (1-omega)*P(nx-1,2) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(nx-2,2) + (1/dy^2)*P(nx-1,3) - R(nx-1,2)); R(nx-1,ny-1) = (rho/dt)*((1/dx)*(u(nx,ny-1) - ufs(nx-2,ny-1)) + (1/dy)*(v(nx-1,ny) - vfs(nx-1,ny-2))); P(nx-1,ny-1) = (1-omega)*P(nx-1,ny-1) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(nx-2,ny-1) + (1/dy^2)*P(nx-1,ny-2) - R(nx-1,ny-1)); P(2,2) = 1.0; E = P - POld; normE = sqrt(sum(sum(E.^2))); error = normE/sqrt(nx*ny); POld = P; if (error >= tol) Iter = Iter + 1; else IterCntr = 1; end end P(1,:) = P(2,:) + (dx/dt)*(u(1,:) - us(1,:)); P(nx,:) = P(nx-1,:) + (dx/dt)*(u(nx,:) - us(nx,:)); P(:,1) = P(:,2) + (dy/dt)*(v(:,1) - vs(:,1)); P(:,ny) = P(:,ny-1) + (dy/dt)*(v(:,ny) - vs(:,ny)); |
|
May 17, 2017, 15:01 |
|
#22 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,783
Rep Power: 71 |
- Why are you prescribing P(2,2)=1??
- Then, you compute E from the difference between p at step k and k-1 but this is not correct...you have to apply the norm on the residual of the equations |
|
May 17, 2017, 15:05 |
P(2,2)
|
#23 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
My professor in my first CFD course said to apply P(2,2) in the corner as a degree of freedom. This has to do with the problem of a singular solution from the PPE. In effect it does not have to be P(2,2) = 1 it could be anything, such that P(2,2) = 5000. Ideally one would specify the average pressure, but this would be computationally expensive. Thank you for spotting the problem in my error.
|
|
May 17, 2017, 15:32 |
|
#24 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,783
Rep Power: 71 |
Quote:
Well, I strongly suggest to cancel such instruction...While it is correct that you have a solution apart from an arbitrary constant value (actually a function of time), fixing a value is not a necessary condition for the convergence. If you write correctly the BC.s in your code, the compatibility condition is fulfilled so that the constant value will be automatically fixed by your SOR iteration. This is also a check that your BC.s are correctly implemented. Furthermore, fixing an arbitrary value would degrade the convergence velocity. |
||
May 18, 2017, 00:01 |
|
#25 | |
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,278
Rep Power: 34 |
Quote:
|
||
May 18, 2017, 13:26 |
Singular solution
|
#26 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Hi,
Thank you for the informative replies. I have implemented my Neumann BC for P(2,2). For my first iteration I get a singularity but after the first iteration I don't get any singularities. Is that supposed to happen? I am under the impression SOR is fixing the singularity after the first iteration. |
|
May 18, 2017, 16:16 |
|
#27 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,783
Rep Power: 71 |
Quote:
|
||
May 19, 2017, 14:19 |
First iteration
|
#28 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
My first iteration gives me the message
Code:
Warning: Contour not rendered for constant ZData > In contourf (line 60) |
|
May 19, 2017, 14:23 |
|
#29 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,783
Rep Power: 71 |
Quote:
Maybe you start you initial guess vector with zeros? However, if you get convergence to the correct solution that's not a problem but (in order to accelerate convergence) consider to save the solution of the previous time step and using it as guess vector for the next time step. |
||
May 22, 2017, 11:42 |
Pressure
|
#30 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
You are correct I am initializing with zeros, with that said I get correct results. One last question in this thread:
I have read both your papers fairly thoroughly and on the boundaries i=1,nx and j=1,ny I am using the non-homogeneous Neumann BC of the following form: This is obviously a better choice than the homogeneous Neumann BC as it satisfies the compatibility condition as discussed in your paper. With that said I do the following Code:
P(1,:) = P(2,:) + (dx/dt)*(us(1,:) - u(1,:)); P(nx,:) = P(nx-1,:) - (dx/dt)*(us(nx,:) - u(nx,:)); P(:,1) = P(:,2) + (dy/dt)*(vs(:,1) - v(:,1)); P(:,ny) = P(:,ny-1) - (dy/dtt)*(vs(:,ny) - v(:,ny)); |
|
May 22, 2017, 11:47 |
|
#31 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,783
Rep Power: 71 |
Quote:
Ok, this way you complete the pressure matrix entries after convergence in the SOR is reached, right? I think that this boundary values will be used also in the correction step. Just check the instruction I highlighted in bold |
|
May 22, 2017, 11:54 |
Pressure BC
|
#32 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Hi Dr. Denaro,
That was a mistake it should be "dt." Furthermore you are correct in saying that it is used in the correction step. Thank you for the immense help you have provided! . |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Centrifugal fan | j0hnny | CFX | 13 | October 1, 2019 13:55 |
Radiation in semi-transparent media with surface-to-surface model? | mpeppels | CFX | 11 | August 22, 2019 07:30 |
My radial inflow turbine | Abo Anas | CFX | 27 | May 11, 2018 01:44 |
Wrong multiphase flow at rotating interface | Sanyo | CFX | 14 | February 7, 2017 17:19 |
Velocity vector in impeller passage | ngoc_tran_bao | CFX | 24 | May 3, 2016 21:16 |