Region couple patch - 1D example
 

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

My results in 1D (stefan1D.tar.gz) track those of BenK in the thread http://www.cfd-online.com/Forums/ope...s-1-6-ext.html (and some others). In particular the patches should be detached to calculate snGrad() and you need to use harmonic interpolation when there is a jump in permeability or the pressure field near the wall is incorrect (nothing new there but just to reiterate).

OpenFOAM assumes the permeability field (or diffusivity) is smoothly varying but there are many practical situations when there is a jump discontinuity at the interface. As I understand it, regionCouplePatch calculates a boundary value of the permeability and then matches the flux from both sides at the boundary. To accomplish this, the pressure on the interface must be incorrect by an amount that depends on the magnitude of the discontinuity in perm and there is also a mismatch in snGrad() on each boundary (both of which can be controlled by increasing the grid resolution). This is exactly what Ben describes.

My question is why not implement the jump discontinuity - in other words match K*snGrad(p) on each side of the boundary with the values of K determined solely from the fields in each domain rather than by coupling them across the patch. I tried two things neither of which worked.
(i) Not updating the boundary conditions on K. This did not help - the K fields remain unchanged but the solver still calculates an average value on the boundary as described by Hrv in the above thread (#20)
(ii) Changing the boundary condition on K from regionCouple to zeroGradient. The solver then failed to converge so that must be wrong.
Perhaps the outcome could depend on whether or not there is a call to K.correctBoundaryConditions().

Before we try to dig into regionCouplePatch, has anyone tried something similar. I see a lot of posts on this topic; Ben seems to have taken it the farthest from what I have found, but still without resolution.

In the 1D case the solution is sufficiently accurate, but I have issues in 2D (see next post) which may be related to this problem.

Regards

Tony

2D Stefan problem
 

My more pressing problem comes with a 2D case. Here the interface between the regions has a cosine profile x=1+A*cos(2*pi*y) with -0.5<y<0.5. For a small amplitude (A << 1) you can solve the coupled Laplace problem by perturbing about the flat interface. So p = B*x + C*exp(+/- u*x)*cos(u*y) with u = 2*pi. I wanted to check that the concentration gradient at the interface (which controls the growth rate) agrees with perturbation theory for small enough amplitudes.
Results comparing OpenFOAM calculations of pressure and concentration with perturbation theory are in plots.tar.gz

With an amplitude A=0.015 the agreement is quite good (*A=0.15*.png). It can be seen that the pressure is converging (albeit a bit slowly - see post #1) to the perturbation solution while the concentration looks to remain a little off.

Since the perturbation analysis is linear we need a double limit of resolution and amplitude so I next tried with A = 0.01, and here is my problem. The variation in pressure just disappears (or more accurately drops by more than an order of magnitude). The agreement does not improve at all with grid resolution (p_A=0.10_N=*.png). Also the change is very sharp. A=0.014 behaves as all the larger amplitudes while A=0.013 behaves like all the smaller amplitudes. I do not understand why there is such an abrupt change in behavior.

The attached case file (stefan2D.tar.gz) should run the case with A=0.015 and N=128 just by typing Allrun. The solver is the same as in #1. You can make other input files using the functions in scripts/tools.py: genPoints() to make a set of boundary points with a cosine profile, newpoints() to make a new blockMeshDict for each region, and linear() to read the data and make the plots.

I would much appreciate any thoughts or comments - I am stuck!

Regards

Tony