[WatchFirefly]

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! (

[bjohn]

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

[WatchFirefly]

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.

[bjohn]

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;

[WatchFirefly]

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.

[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

[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.

[WatchFirefly]

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,

[bjohn]

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

[WatchFirefly]

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.

[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.

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!

[bjohn]

>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.

[WatchFirefly]

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.

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.