|
[Sponsors] |
Help with Simple 2D Incomp. Inviscid Unsteady FDM Solver |
![]() |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
![]() |
![]() |
#1 |
New Member
David Hilbert
Join Date: May 2024
Posts: 1
Rep Power: 0 ![]() |
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: 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); } } } Last edited by David Hilbert; May 28, 2024 at 05:29. Reason: Images didn't work the first time |
|
![]() |
![]() |
![]() |
![]() |
#2 |
New Member
Hardware
Join Date: May 2024
Posts: 7
Rep Power: 2 ![]() |
1.Smaller time step
2.No-normal-flux pressure BCs 3.Double precision calculations Visit Frelan Hardware for burlington handles Last edited by michealkors; June 4, 2024 at 04:47. |
|
![]() |
![]() |
![]() |
![]() |
#3 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
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. |
|
![]() |
![]() |
![]() |
Tags |
beginner, fdm, projection method |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
Time dependent solver of unsteady Navier Stokes | Mehrez | COMSOL | 10 | March 20, 2024 03:57 |
Divergence detected in AMG solver | Rahmani_Zakaria | Fluent Multiphase | 4 | January 11, 2023 15:09 |
Unsteady Restart Divergence | pro_ | SU2 | 6 | May 20, 2020 15:17 |
[ANSYS Meshing] Help with element size | sandri_92 | ANSYS Meshing & Geometry | 14 | November 14, 2018 07:54 |
comments on FDM, FEM, FVM, SM, SEM, DSEM, BEM | kenn | Main CFD Forum | 2 | July 18, 2004 18:28 |