[David Hilbert]

Hello, I am new here. I am currently taking a course on CFD, so I was interested in putting the theory into practice and making the simplest possible solver for a 2D problem unsteady problem. In a basic way, it sometimes works, but I am definitely having problems with stability and just the general result not being what it should be, so, after trying many things to fix it, I concluded my best option is to ask here.




My approach was to use Chorin's Projection Method (https://en.wikipedia.org/wiki/Projec...luid_dynamics)) along with FDM, which seems to theoretically deal with the nonlinearity and the implicit pressure term without having to do anything fancy.

I used Unity to do the visualizing and such, because I use it for all sorts of things. The language is c#, and my code is in attachments, FDM_Incomp_Inviscid_2D_chorinsProj.txt. (CFD Online does not accept .cs attachment)

So, there is a lot of unimportant stuff that cannot be the problem. The actual calculations are from line 188-319, especially the time loop (I also put it at the end in a code block). I will try to explain what I am doing.

Initial Conditions:
For my initial conditions, I have tried various things: I have tried a small region of non-zero velocity in the center of the grid, or velocity completely 0, and I have tried to make pressure either zero or set the analytical hydrostatic pressure. Neither improve the result, so for simplicity I go with everything = 0, and instead introduce some unsteadiness through the body force (gravity) and time-dependent boundary conditions.

Within the time loop:
Step 1. Calculate intermediate velocity U_star, V_star ("prediction step")
From Chorin's Projection method, the first step is to calculate the intermediate velocity using only values from the previous time step, and neglecting the pressure term.

https://ibb.co/mypvx3Y

I use the second order finite differences for the spatial discretization in lines 247, 248, and I also introduce the same boundary conditions for U_star and V_star as I would like for U^{n+1} and V^{n+1}.


Step 2. Solve Pressure Poisson Equation
To get the Pressure Poisson Equation, first you use the "projection step".

https://ibb.co/Z21zMhk

As I understand it, the reason this more or less works is because if you plug it back into the previous equation, the one relating u_star and u^n, you get the actual momentum conservation equation with forward Euler time step, but also using p^{n+1} rather than p^n. The fact that this method exists tells me that this is valid and can converge.

Ok so div of that equation and you get the Pressure Poisson Equation

https://ibb.co/k9Kp9Lg
which I solve iteratively in lines 252-293 using the Jacobi method because it is the simplest way. I also impose the boundary conditions of grad p in normal direction of the boundary = 0.




Step 3. Get u^{n+1} from the projection step

https://ibb.co/Z21zMhk
Using this equation again, this time because we have u_star and p^{n+1}, we can solve for u^{n+1}. Note that in my code, in the time loop, n is shifted back by one compared to the wikipedia formulas.


Now, here are some of my simulation tests (CFD Online does not accept .mp4 attachment):
Cavity Flow:
https://youtu.be/no84Cnh6mbw
https://youtu.be/ZfXuIJZIBik
Pipe Flow (yes the BC is simply uniform on both sides..not realistic):
https://youtu.be/3f0WMlEJCiU
Just gravity:
https://youtu.be/8B9JJRvLw-s


Specific information in the video descriptions.



I would like to note that one of the biggest problems other than the eventual instability of the velocity (not seen in some of the videos because it didn't run long enough) is that the pressure decreases everywhere as time goes on, which is not physical. My example of just gravity simulation only shows the gradient well because I did a very short time interval. But for longer simulations the gradient is pretty much invisible compared to the overall change (decrease) of pressure. (For visualizing the pressure I take the fraction (p-pMin)/(pMax/pMin), where the extrema are from any point at any time)


Another specific obervation is that, in the just gravity simulation, the velocity doesn't actually converge to 0--except for every second row for some reason..?


Anyways, thanks in advance for any help you can give. I hope I am doing everything right with the posting, sometimes forums have specific rules and such idk.





For easy reference, here is the time loop portion of the code:

// // time loop
        double dx2 = Pow(dx, 2); // calculated only once for efficiency
        double dy2 = Pow(dy, 2); //
        for (int n = 1; n < nCount; n++)
        {
            // // 1. intermediate velocity
            double[,] U_star = new double[xCount, yCount];
            double[,] V_star = new double[xCount, yCount];

            //BCs
            for (int i = 1; i < xCount - 1; i++) // top and bottom
            {
                U_star[i, yCount - 1] = U_onUpperBoundary(xBounds[0] + i * dx, n * dt);
                V_star[i, yCount - 1] = 0;
                U_star[i, 0] = 0;
                V_star[i, 0] = 0;
            }
            for (int j = 1; j < yCount - 1; j++) // right and left
            {
                U_star[xCount - 1, j] = 0;
                V_star[xCount - 1, j] = 0;
                U_star[0, j] = 0;
                V_star[0, j] = 0;
            }
            for (int i =1; i < xCount-1; i++)
            {
                for (int j =1; j < yCount-1; j++)
                {
                     U_star[i, j] = U[i, j, n - 1] + dt * (-(U[i, j, n - 1] * (U[i + 1, j, n  - 1] - U[i - 1, j, n - 1]) / (2 * dx) + V[i, j, n - 1] * (U[i, j + 1, n  - 1] - U[i, j - 1, n - 1]) / (2 * dy)) + bodyForce[0]);
                     V_star[i, j] = V[i, j, n - 1] + dt * (-(U[i, j, n - 1] * (V[i + 1, j, n  - 1] - V[i - 1, j, n - 1]) / (2 * dx) + V[i, j, n - 1] * (V[i, j + 1, n  - 1] - V[i, j - 1, n - 1]) / (2 * dy)) + bodyForce[1]);
                }
            }

            // // 2. iterate poisson eq
            // initial pressure guess = pressure from time step n-1
            for (int i = 0; i < xCount; i++)
            {
                for (int j = 0; j < yCount; j++)
                {
                    P[i, j, n] = P[i, j, n - 1];
                }
            }
            // iterate
            double[,] pNew = new double[xCount, yCount];
            for (int iter = 0; iter < iterMax; iter++)
            {
                for (int i = 1; i < xCount - 1; i++)
                {
                    for (int j = 1; j < yCount - 1; j++)
                    {
                         pNew[i, j] = 1 / (2* (dx2 + dy2)) * ((P[i + 1, j, n] + P[i - 1, j, n]) *  dy2 + (P[i, j + 1, n] + P[i, j - 1, n]) * dx2 -
                             dx2 * dy2 * density / dt * ((U_star[i + 1, j] - U_star[i - 1, j]) / (2 *  dx) + (V_star[i, j + 1] - V_star[i, j - 1]) / (2 * dy)));
                    }
                }

                // pressure BCs
                for (int i = 1; i < xCount - 1; i++) // top and bottom impermeable
                {
                    pNew[i, 0] = pNew[i, 1];
                    pNew[i, yCount - 1] = pNew[i, yCount - 2];
                }
                for (int j = 1; j < yCount - 1; j++) // right and left impermeable
                {
                    pNew[0, j] = pNew[1, j];
                    pNew[xCount - 1, j] = pNew[xCount - 2, j];
                }

                for (int i = 0; i < xCount; i++)
                {
                    for (int j = 0; j < yCount; j++)
                    {
                        P[i, j, n] = pNew[i, j];
                    }
                }
            }

            // // 3. get u^{n+1}
            //BCs
            for (int i = 1; i < xCount - 1; i++) // top and bottom
            {
                U[i, yCount - 1, n] = U_onUpperBoundary(xBounds[0] + i * dx, n * dt);
                V[i, yCount - 1, n] = 0;
                U[i, 0, n] = 0;
                V[i, 0, n] = 0;
            }
            for (int j = 1; j < yCount - 1; j++) // right and left
            {
                U[xCount - 1, j, n] = 0;
                V[xCount - 1, j, n] = 0;
                U[0, j, n] = 0;
                V[0, j, n] = 0;
            }
            for (int i = 1; i < xCount-1; i++)
            {
                for (int j = 1; j < yCount-1; j++)
                {
                    U[i, j, n] = U_star[i, j] - dt / density * (P[i + 1, j, n] - P[i - 1, j, n]) / (2 * dx);
                    V[i, j, n] = V_star[i, j] - dt / density * (P[i, j + 1, n] - P[i, j - 1, n]) / (2 * dy);
                }
            }
        }

[michealkors]

1.Smaller time step
2.No-normal-flux pressure BCs
3.Double precision calculations

Visit Frelan Hardware for burlington handles

[FMDenaro]

I'm not going to read your code, but:


1) BCs for the intermediate velocity are required only in case of implicit time-integration, not for explicit methods.
2) Use a staggered grid for the variables. This is the basis of the Exact Projection Method.
3) The pressure equation must be closed by proper Neumann BCs that satisfy the compatibility condition to ensure the existence of one solution.
4) At the end of the first time step check the divergence-free condition in each cell, it must be zero at the precision of the iterative solver.