Problem with validation of regionCouple solver
 

I am working on a solver for Stefan problems - i.e. a sharp interface (typically moving) separating two different regions. I started with conjugateHeatFoam from foam-extend-3.0. I solve a Laplace equation for the pressure in two regions with different permeability and then solve for a concentration field (convection-diffusion) in the inner region. Solver is attached (stefanFoam.tar.gz).

In the test problem the interface between the regions has a cosine profile x=1+A*cos(2*pi*y) with -0.5<y<0.5, with A constant. For clarity I choose both perms to be unity so that on the boundary between the domains the pressure should be simply -x if the pressure at the inlet is zero and the pressure gradient at the outlet is -1. The plots show the pressure from OpenFOAM (blue) vs the analytic solution (red). If the amplitude A is 0.02 (or larger) the solution is correct but for A=0.01 the variation in pressure is much too small (by more than an order of magnitude). In fact there is a more or less discontinuous change in behavior; A=0.014 is fine and A=0.013 is not. I have played with grid resolution, domain size, fvSchemes and fvSolution but to no avail.

The illustrated problems can be run from the attached tarballs. The python functions in scripts/tools/py can be used to make new blockMeshDict files if needed.

Does anyone have any suggestions?

Regards

Tony