2D Finite Volume Method  Please help
Hi Everyone,
I am trying to solve a 2D shocktube problem in a 2D, cellcentered, First order finite volume framework. I coded Roe's scheme using the sample code given at I do like cfd site, but after I run the code nothing seems to happen to density and pressure, even though I let the code run for really long time. I don't see any initial structure (shock, contact, expansion) at the interface of left and right sections of the grid. Whle for u and v, there're some really strong disturbances only in some cells. I have a feeling that either I have not set up my normal vectors correctly or that my choice of left and right states at a given interface is not correct. * In the first order cellcentered FVM, we define values at the center of the cell, and use the same values at the cellinterface. In secondorder, we will reconstruct the variation of primary variables to obtain their interface values from cellcenters. Is that correct? * Now, in a 2D FVM method, we solve Riemann Problem at every interface, right? Consider a left face of cell (i,j). Now the When I am solving a Riemann problem on this left interface, is the following correct? Cell : (i,j) Interface: Left Left state for above interface: Left Interface Value of cell (i,j) Right state: Right Interface value of Cell (i1,j) Similarly, for the same cell, but Interface: Bottom Left state: Bottom Interface Value of cell (i,j) Right State: Top Interface Value of cell (i,j1) Is above correct? Thank you very much. Please let me know if you need further info. PS: I am actually kinda very frustrated b/c, I have am stuck at 2D FVM, what a shame! :(( 
I think it's correct that you say.
Some things to check: 1) Add up face normals for each cell. The result should be zero. 2) Set up a uniform flow over the domain and run the code. The solution should not change. 3) Send the left and the right states of the initial discontinuity into the Roe flux subroutine to see if the flux has zero values in the 1st and the 4th components (density and pressure). If it does, the flux subroutine is not correct. It is a good idea to focus on one cell (at the initial discontinuity) and investigate why the density and pressure do not change. Good luck. bjohn 
Hi Bjohn,
Thanks for a quick reply. 1) I checked face normals. I added up Δx, and Δy for all four faces of a given cell, and that sum came out to be zero for all the cells of this 20 x 20 rectangular grid. I have defined Δx and Δy as follows: (going anticlockwise) //Right face deltaxR[i][j] = x[i+1][j+1]  x[i+1][j]; deltayR[i][j] = y[i+1][j+1]  y[i+1][j]; //Top deltaxT[i][j] = x[i][j+1]  x[i+1][j+1]; deltayT[i][j] = y[i][j+1]  y[i+1][j+1]; //Left deltaxL[i][j] = x[i][j]  x[i][j+1]; deltayL[i][j] = y[i][j]  y[i][j+1]; //Bottom deltaxB[i][j] = x[i+1][j]  x[i][j]; deltayB[i][j] = y[i+1][j]  y[i][j]; I also added up the normalized Δx, and Δy. Their sum too came out to be zero. Normalized Δx = Δx/Length of the respective face. 2) I set up the uniform flow as you suggested, and surprisingly, the value of v (ycomponent of velocity), had some value instead of 0.0. Velocity components:
In Roeflux subroutine, Left flux is defined as follows: f1Left = rhoLeft*qLeft; f2Left = rhoLeft*qLeft*uLeft + pLeft*ny; f3Left = rhoLeft*qLeft*vLeft + pLeft*nx; f4Left = rhoLeft*qLeft*HLeft; Right flux f1Right = rhoRight*qRight; f2Right = rhoRight*qRight*uRight + pRight*ny; f3Right = rhoRight*qRight*vRight  pRight*nx; (this was + initially.) f4Right = rhoRight*qRight*HRight; Only when I made the + sign  in f3Right, I got the v to be 0.0 (as set up initially). This is really weird because f3right and f3left has different formulas!!!! I need to check my roeflux routine again. Bjohn, please stay with me on this one. I want to really finish this project. It's been pending for a long time. Thanks, and please let me know your comments. Once, I have checked my roeflux routine, I will post a message. 
1) Your (deltax,deltay) is tangent to the face. If you define deltax and deltay as you describe, normal will be (deltay, deltax).
2) The vector multiplied by pressure should be the normal: f2Left = rhoLeft*qLeft*uLeft + pLeft*nx; f3Left = rhoLeft*qLeft*vLeft + pLeft*ny; Right flux f2Right = rhoRight*qRight*uRight + pRight*nx; f3Right = rhoRight*qRight*vRight + pRight*ny; 
Bjohn,

>qR[i][j] = u[i][j]*deltayR[i][j]v[i][j]*deltaxR[i][j]; //Normal and
>rR[i][j] = u[i][j]*deltaxR[i][j]+v[i][j]*deltayR[i][j]; //tangent. These quantities must be based on the unit vectors: ds = sqrt(deltaxR[i][j]^2 + deltayR[i][j]^2) qR[i][j] = ( u[i][j]*deltayR[i][j]v[i][j]*deltaxR[i][j] )/ds rR[i][j] = ( u[i][j]*deltaxR[i][j]+v[i][j]*deltayR[i][j] )/ds >When I ran the code for one iteration, I noticed the the sum of third component of >normal flux is NOT zero, and the v = 0.0004!! It should be zero! Nothing is moving. This is not good. >Then, I changed the sign in front of pRight from + to  in f3Right, and value of v >became zero as expected. Is this correct? No, it is not. The pressure term should be based on the normal: +pressure*nx +pressure*ny not +pressure*ny +pressure*nx nor +pressure*ny pressure*nx Normal vector should be the unit vector inside the flux subroutine (unless it is normalized in the subroutine). Typically, you multiply the Roe flux by the magnitude of the face (ds) at the end. bjohn 
>ρ_R = ρ_L = 1.0, p_L = p_R = 1.0, u = v = 0.0 through out.
In this case, you have f1Left = 0 f2Left = nx f3Left = ny f4Left = 0 Right flux f1Right = 0 f2Right = nx f3Right = ny f4Right = 0 and the dissipation term should be zero since Delta_U = 0. So, the Roe flux becomes f1 = 0 f2 = nx f3 = ny f4 = 0 If your normals add up to zero within a cell, the residuals will be zero, and thus the solution should not change at all. 
Okay. Gotcha.
I will check my code for again for everything that you mentioned. My code is based on a sample code given at the site of I do like cfd book. http://www.ossanworld.com/cfdbooks/c..._fluxes_v2.f90 My code is in C, and the funny thing is even tho I have the above code to refer to my code still does not work. Must make it work. How about the initial values for a test case? and Is my expectation to see that expansioncontactshock structure in 2D correct? thanks, 
The sample code looks good to me.
Maybe, errors have been introduced when you convert it into C code... Looks like you're trying to run Sod's shock tube problem. There is a plot in the Wiki page: http://en.wikipedia.org/wiki/Sod_shock_tube bjohn 
Yep. Thanks for the link. So I should basically expect the same structure, but in 2D, so plannar wave, basically.
I will update you soon...tonight.. >>Edit>> Also, wanted to mention that I started writing code before I checked it out on cfdbooks.com (as a last resort), and my implementation is a bit different. Though cfdbooks.com's code is kinda cleaner than mine! Post my findings later.. man, I gotta finish this 2D and then move on to writing 3D codes. 
bJohn,
One of the major error I found is that my Δx and Δy are just the coordinate differences, and not absolute Δx values. So, my nx, and ny are directly Δy, and Δx without negative sign for Δy. Once I corrected that I started get zero flux for the uniform condition flow I posted earlier. I am still checking my code, by setting up u_L=1.0 and u_R=1.0 i.e. two opposing flows at the discontinuity. With pressure being the same, this should be a standing flow. Though this does not sound physically possible. Anyways, I will post my findings later. Let me know if you have any suggestions. Thanks, man! 
>u_L=1.0 and u_R=1.0 i.e. two opposing flows at the discontinuity
This is probably not a steady shock satisfying the RankineHugoniot relation. So, waves will come out. Sod' shock tube problem mentioned above is a good test case. It can be used in various directions (leftright, rightleft, topbottom, bottomtop, or even oblique directions) to see if the code is able to produce the same wave structure (expansioncontactshock) independently of the direction. 
All times are GMT 4. The time now is 08:00. 