Register Blogs Members List Search Today's Posts Mark Forums Read

 September 4, 2010, 22:10 2D Finite Volume Method - Please help #1 New Member   Watch Join Date: Sep 2010 Posts: 10 Rep Power: 8 Hi Everyone, I am trying to solve a 2D shocktube problem in a 2D, cell-centered, 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 cell-centered FVM, we define values at the center of the cell, and use the same values at the cell-interface. In second-order, we will reconstruct the variation of primary variables to obtain their interface values from cell-centers. 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 (i-1,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,j-1) 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: 8 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: 8 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 anti-clockwise) //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 (y-component of velocity), had some value instead of 0.0. Velocity components: u = x-compo v = y-compo q = normal to face velo. r = tangential to face velo. In Roe-flux 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 roe-flux 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 roe-flux routine, I will post a message.

 September 6, 2010, 06:26 #4 New Member   Join Date: Jun 2010 Posts: 16 Rep Power: 8 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; tienthanh_bk2901 likes this.

 September 8, 2010, 00:10 #5 New Member   Watch Join Date: Sep 2010 Posts: 10 Rep Power: 8 Bjohn, My normal and tagents are defined the way you mentioned. e.g. 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. I set up the initial condition for the uniform flow as follows: ρ_R = ρ_L = 1.0, p_L = p_R = 1.0, u = v = 0.0 through out. 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. Then, I changed the sign in front of pRight from + to - in f3Right, and value of v became zero as expected. Is this correct? After this, I ran with ic as given in item-8, no joy. I just see standing shock in the density plot, and corners are rounded that's it. Also, do you have results of 2D shock tube on 20 x 20 rectangular grid, with domain of 8 x 8. I just want to see the results, so that I know what to expect. I am currently expecting to see same expansion, contact, shock structure as we see in 1D shock tube problem. What initial conditions should I use? I was using ρ_R = 0.125, ρ_L = 1.0, p_L = 0.1, p_R = 1.0, u = v = 0.0 Thanks, man. I appreciate.

 September 8, 2010, 05:02 #6 New Member   Join Date: Jun 2010 Posts: 16 Rep Power: 8 >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 tienthanh_bk2901 likes this. 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: 8 >ρ_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. tienthanh_bk2901 likes this.

 September 8, 2010, 10:20 #8 New Member   Watch Join Date: Sep 2010 Posts: 10 Rep Power: 8 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 expansion-contact-shock structure in 2D correct? thanks,

 September 8, 2010, 10:57 #9 New Member   Join Date: Jun 2010 Posts: 16 Rep Power: 8 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: 8 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: 8 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: 8 >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 Rankine-Hugoniot relation. So, waves will come out. Sod' shock tube problem mentioned above is a good test case. It can be used in various directions (left-right, right-left, top-bottom, bottom-top, or even oblique directions) to see if the code is able to produce the same wave structure (expansion-contact-shock) independently of the direction.

September 12, 2010, 20:58
#13
New Member

Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 8
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 anti-clockwise.

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.x-momentum 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. x-momentum 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:
 Originally Posted by WatchFirefly 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.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post falopsy Main CFD Forum 44 January 15, 2014 03:17 Scott2 Main CFD Forum 1 August 16, 2010 20:24 littlelz Main CFD Forum 5 June 21, 2003 12:50 Aditya Vaze Main CFD Forum 1 January 19, 2000 14:55 Anthony Main CFD Forum 4 February 2, 1999 20:24

All times are GMT -4. The time now is 08:31.