CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

2D Finite Volume Method - Please help

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

Like Tree3Likes
  • 1 Post By bjohn
  • 1 Post By bjohn
  • 1 Post By bjohn

Reply
 
LinkBack Thread Tools Display Modes
Old   September 4, 2010, 22:10
Unhappy 2D Finite Volume Method - Please help
  #1
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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! (
WatchFirefly is offline   Reply With Quote

Old   September 4, 2010, 23:25
Default
  #2
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
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
bjohn is offline   Reply With Quote

Old   September 5, 2010, 22:41
Default
  #3
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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.
WatchFirefly is offline   Reply With Quote

Old   September 6, 2010, 06:26
Default
  #4
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
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.
bjohn is offline   Reply With Quote

Old   September 8, 2010, 00:10
Default
  #5
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
Bjohn,

  1. 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.
  1. I set up the initial condition for the uniform flow as follows:
  2. ρ_R = ρ_L = 1.0, p_L = p_R = 1.0, u = v = 0.0 through out.
  3. 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.
  4. Then, I changed the sign in front of pRight from + to - in f3Right, and value of v became zero as expected. Is this correct?
  5. 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.
  6. 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.
  7. 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.
WatchFirefly is offline   Reply With Quote

Old   September 8, 2010, 05:02
Default
  #6
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
>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.
bjohn is offline   Reply With Quote

Old   September 8, 2010, 05:50
Default
  #7
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
>ρ_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.
bjohn is offline   Reply With Quote

Old   September 8, 2010, 10:20
Default
  #8
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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,
WatchFirefly is offline   Reply With Quote

Old   September 8, 2010, 10:57
Default
  #9
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
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
bjohn is offline   Reply With Quote

Old   September 8, 2010, 11:06
Default
  #10
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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.
WatchFirefly is offline   Reply With Quote

Old   September 10, 2010, 22:43
Default
  #11
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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!
WatchFirefly is offline   Reply With Quote

Old   September 10, 2010, 23:08
Default
  #12
New Member
 
Join Date: Jun 2010
Posts: 16
Rep Power: 6
bjohn is on a distinguished road
>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.
bjohn is offline   Reply With Quote

Old   September 12, 2010, 20:58
Arrow
  #13
New Member
 
Watch
Join Date: Sep 2010
Posts: 10
Rep Power: 5
WatchFirefly is on a distinguished road
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 (n_x,ny), then

n_x = \Delta y
n_y = -\Delta x

where,

\Delta x = x_{i+1,j+1}-x_{i+1,j} No need to use absolute values.
\Delta y = y_{i+1,j+1}-y_{i+1,j}

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:
u_L = u_R = v_L =v_R=0.0
\rho_L = 1.0
\rho_R = 0.125
p_L = 10.0
p_R= 1.0

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. f_1^*, f_2^*,f_3^*, f_4^* for both Left and Right state.

At the end of the first iteration,
\sum f_1^* > 0 i.e. mass flux was positive
\bf{\sum f_2^* < 0} 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.
{u2_{ij}}^{n+1}   = {u2_{ij}}^{n} - \frac{dt}{A_{ij}} \left( Fn2R_{ij}+Fn2T_{ij}+Fn2L_{ij}+Fn2B_{ij}\right);
where,
Fn2 = f_2^*

\sum f_3^* = 0, as expected.
\sum f_4^* > 0 as expected.

PS: Because of your help, I am a lot more motivated to finish this project, so thanks.





Quote:
Originally Posted by WatchFirefly View Post
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.
WatchFirefly is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Finite Volume Method For Cylindrical Coordinates falopsy Main CFD Forum 44 January 15, 2014 02: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 13:55
Data Structure for the unstructured finite volume method Anthony Main CFD Forum 4 February 2, 1999 19:24


All times are GMT -4. The time now is 05:07.