[selig5576]

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.

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

I have run a simulation with 100^2 points until T = 12.0 and my results seem to be very good. I have attached photos to illustrate.

[FMDenaro]

- 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

[selig5576]

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.

[FMDenaro]

Originally Posted by selig5576

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.

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.

[arjun]

Originally Posted by FMDenaro

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.

Fixing also has other numerical issues of convergence related.

[selig5576]

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.

[FMDenaro]

Originally Posted by selig5576

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.

What do you mean for "singularity "?? You get NaN??

[selig5576]

My first iteration gives me the message

Warning: Contour not rendered for constant ZData 
> In contourf (line 60)

The first iteration gives a P full of zeros, not NaNs, so I don't think its too much of an issue. If I have P(2,2) = 1.0 instead of the corner Neumann BC I get a similar result, so I don't think its too much of an issue?

[FMDenaro]

Originally Posted by selig5576

My first iteration gives me the message

Warning: Contour not rendered for constant ZData 
> In contourf (line 60)

The first iteration gives a P full of zeros, not NaNs, so I don't think its too much of an issue. If I have P(2,2) = 1.0 instead of the corner Neumann BC I get a similar result, so I don't think its too much of an issue?

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.

[selig5576]

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

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

In your paper you do a similar discretization process except your h is non assuming uniformity like I am (for now.) Note, I apply these BCs after my iteration process.

[FMDenaro]

Originally Posted by selig5576

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

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

In your paper you do a similar discretization process except your h is non assuming uniformity like I am (for now.) Note, I apply these BCs after my iteration process.

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

[selig5576]

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