[kupiainen]

Is there a reflecting boundary condition in OpenFoam?

[henry]

fixedValue, zeroGradient, symmetryPlane and slip boundary conditions are all reflecting.

[kupiainen]

I do not mean reflecting as in not non-reflecting, but as in reflecting everything i.e. (F,n) = 0,
where F is the flux and n is the normal.

[henry]

Isn't that a symmetryPlane?

[kupiainen]

Well I don't know what a symmetryPlane is mathematically, but it seems that it doesn't reflect everything. I run the Woodward, Colella (JCP 54, 1984) forwardfacing step testcase, and the top and bottom boundaries should be set to reflecting. I run on a very fine mesh O(1e5 cells) and I expect to see a contact discontinuity near the top boundary, but it is not there. I believe that the contact discontinuity is due to the reflecting b.c. Right now I am using symmetryPlanes as boundary conditions.

BTW the OpenFoam manual could be little more specific about what the b.c. really mean.

[henry]

I run with symmetryPlanes and the shocks reflect perfectly well.

The OpenFOAM manual will be more specific about lots of things as soon as we obtain funding in the form of donations/sponsorship to cover the cost of writing it. So far we have received ~$300 in donations and none in sponsorship so we are not in a position to continue writing documentation.

[kupiainen]

Do you resolve the contact discontinuity near the top boundary to the left of the bowshock?

[henry]

I haven't studied the results in detail for the old sonic-solver. I have developed a new version of rhopSonicFoam for the 1.2 release which is undergoing extensive testing at the moment as part of a research project, it might be interesting for you to test it on your case.

[kupiainen]

I am looking forward to the 1.2 release.

[kupiainen]

Another thing with this testproblem is that the shock is moving with the wrong speed compared to the results of Colella and Woodward. The shock eventually propagates out of the domain to the left.
I now seem to resolve the contact discontinuity, although it requires a lot of cells.
Any ideas with the incorrect shock speeds?

[henry]

Which code are you running?

[kupiainen]

I do an Euler step in time (can be extended to a TVD Runge-Kutta 2:nd order) and use central differences in space with Jameson artificial dissipation.
The laplacian is discretized using
Gauss linear correction
and divergence terms using
Gauss linear

This is the problably the most classical way of solving this problem, however I am not sure about the reflecting b.c..

[henry]

I was assuming you were running one of the standard OpenFOAM sonic-solvers; I have no idea why your code produces the incorrect shock-speed.

[kupiainen]

The reason can be twofold, either the symmetryPlane is not equivalent to a reflecting b.c. or there is some bug in the div, laplacian etc. operators in forms of conservation. If the operators are not discretized in the correct way, errors of this type might show. For solutions not involving discontinuities these errors will not show. It is my belief that the errors are not due to the operators though, but I will continue searching for answers. So if anybody has some good and constructive ideas I am more than welcome to try them.

[kupiainen]

I have run the rhoSonicFoam code now. In order to run this problem, I had to modify the pressure, since the gas is a gamma-law gas. So
p = (gamma-1)*(rhoE-0.5*rho*magSqr(rhoU/rho));
and the temperature is not needed (and not defined in this setting). The pressure is set to 0 everywhere and velocity is = (3 0 0) everywhere and the density is set to 1.4 everywhere. These are also the inlet b.c. for this problem. At the outflow zeroGradient is used for all components, since it is a supersonic outlet. On all other physical boundaries I use symmetryPlanes, since these are the reflecting b.c. as described in the article by Colella and Woodward.
The code is stable for CFL-number 0.6 at least, but the resuls are puzzeling. At the dimensionless time 4 the results from output from rhoSonicfoam and Colella et. al. are compared. I have refined the grid ad absurdum (up to 500 000 cells) and made grid convergence studies and found that rhoSonicfoam produces extra and new shockwaves, which are not present in the Colella et. al. paper. The main shock seems to move with the right speed, but the smooth part of the solution is very wiggly. The new structures seen in the solution are due to numerical errors and should not be there.

As I understand OpenFoam the variable
phi = linearInterpolate(rho*U & mesh.Sf());
is used to find the upwind direction. When the flow is supersonic the information does not however move only with the flow speed. This is well-known from hyperbolic conservation laws. One would rather need to use a tensorvalued variable for the upwind direction, since it is different for the different components of the e.g. div operator. The basic fact is that one misses the characteristics related to mag(U)-SoundSpeed.
One possible way to circumvent this problem, which is not so elegant in my view, is to use the PISO coupling. Here the pressure iterations make the problem globally coupled. This procedure is natively developed for incompressible flows, where this is natural, since the soundspeed is infinity and no numerical method used for hyperbolic problems works well (characteristics are finite in hyperbolic problems).

Does anybody have another explanation for this? Comments? I am happy to discuss these issues, since I want to solve this problem.

[kupiainen]

Problem solved:
The artificial viscosity was not implemented in conservation form.
Marco