# Region couple patch - 1D example

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

May 22, 2014, 10:28
Region couple patch - 1D example
#1
Member

Join Date: Aug 2013
Posts: 37
Rep Power: 6
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
Attached Files
 stefanFoam.tar.gz (2.8 KB, 5 views) tests.xlsx (5.6 KB, 2 views) stefan1D.tar.gz (16.3 KB, 5 views)

May 22, 2014, 11:24
2D Stefan problem
#2
Member

Join Date: Aug 2013
Posts: 37
Rep Power: 6
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
Attached Files
 stefan2D.tar.gz (22.1 KB, 2 views) plotsC.tar.gz (70.0 KB, 6 views) plotsp-10.tar.gz (64.9 KB, 3 views) plotsp-15.tar.gz (70.3 KB, 4 views) plotsp-1314.tar.gz (61.4 KB, 1 views)

Last edited by tladd; May 22, 2014 at 11:36. Reason: Added files

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post skuznet OpenFOAM Running, Solving & CFD 99 March 16, 2017 06:07 tH3f0rC3 OpenFOAM 7 February 23, 2013 06:37 cfdengineering OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ... 48 January 25, 2013 04:28 matteo_gautero OpenFOAM Running, Solving & CFD 0 February 28, 2008 07:51 michele OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ... 2 July 15, 2005 04:15

All times are GMT -4. The time now is 15:21.