
[Sponsors] 
September 4, 2010, 22:10 
2D Finite Volume Method  Please help

#1 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
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! ( 

September 4, 2010, 23:25 

#2 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
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 

September 5, 2010, 22:41 

#3 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
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. 

September 6, 2010, 06:26 

#4 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
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; 

September 8, 2010, 00:10 

#5 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
Bjohn,


September 8, 2010, 05:02 

#6 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
>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 Last edited by bjohn; September 8, 2010 at 05:59. 

September 8, 2010, 05:50 

#7 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
>ρ_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. 

September 8, 2010, 10:20 

#8 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
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, 

September 8, 2010, 10:57 

#9 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
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 

September 8, 2010, 11:06 

#10 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
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. Last edited by WatchFirefly; September 8, 2010 at 14:55. 

September 10, 2010, 22:43 

#11 
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
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! 

September 10, 2010, 23:08 

#12 
New Member
Join Date: Jun 2010
Posts: 16
Rep Power: 7 
>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. 

September 12, 2010, 20:58 

#13  
New Member
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 6 
John,
What I wrote below about nx and ny is actually not correct. I found that without ve sign in front of Δx, I was not getting uniform flow as I should and later I checked it manually and realized that I, indeed, need to put negative sign in front of Δx's and do not need to use absolute values. so, if a normal vector to the right face is , then where, No need to use absolute values. and so on for other faces, going anticlockwise. After changing my definitions for n_x and n_y, I set up uniform flow, and ran the code for 1 and 10 iterations, got absolutely no change in any values, as expected.  Then, I set up sod's shocktube case as follows: For one of the cells, where I have set up the discontinuity, I was looking at all the primitive variables including normal and tangential velocities, and each component of fluxes of the solutions i.e. for both Left and Right state. At the end of the first iteration, i.e. mass flux was positive i.e.xmomentum flux was negative. Is this surprising? Shouldn't this be positive? Net flux going out should be positive? I should mention that while this sum is negative in the final update formula, the, second component u2 at time n+1, shown below, i.e. xmomentum flux is +ve. where, , as expected. as expected. PS: Because of your help, I am a lot more motivated to finish this project, so thanks. Quote:


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Finite Volume Method For Cylindrical Coordinates  falopsy  Main CFD Forum  44  January 15, 2014 03:17 
Chorin's Projection Method for Finite Volume  Scott2  Main CFD Forum  1  August 16, 2010 20:24 
ALE in finite volume method  littlelz  Main CFD Forum  5  June 21, 2003 12:50 
finite volume method for CFD analysis of 2D blunt body  Aditya Vaze  Main CFD Forum  1  January 19, 2000 14:55 
Data Structure for the unstructured finite volume method  Anthony  Main CFD Forum  4  February 2, 1999 20:24 