CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Running, Solving & CFD

Region couple patch - 1D example

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

Like Tree3Likes
  • 1 Post By tladd
  • 2 Post By tladd

LinkBack Thread Tools Display Modes
Old   May 22, 2014, 10:28
Default Region couple patch - 1D example
Tony Ladd
Join Date: Aug 2013
Posts: 37
Rep Power: 6
tladd is on a distinguished road
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 (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.


Attached Files
File Type: gz stefanFoam.tar.gz (2.8 KB, 5 views)
File Type: xlsx tests.xlsx (5.6 KB, 2 views)
File Type: gz stefan1D.tar.gz (16.3 KB, 5 views)
mm.abdollahzadeh likes this.
tladd is offline   Reply With Quote

Old   May 22, 2014, 11:24
Default 2D Stefan problem
Tony Ladd
Join Date: Aug 2013
Posts: 37
Rep Power: 6
tladd is on a distinguished road
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/ 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!


Attached Files
File Type: gz stefan2D.tar.gz (22.1 KB, 2 views)
File Type: gz plotsC.tar.gz (70.0 KB, 6 views)
File Type: gz plotsp-10.tar.gz (64.9 KB, 3 views)
File Type: gz plotsp-15.tar.gz (70.3 KB, 4 views)
File Type: gz plotsp-1314.tar.gz (61.4 KB, 1 views)

Last edited by tladd; May 22, 2014 at 11:36. Reason: Added files
tladd is offline   Reply With Quote


Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
conjugate heat transfer in OpenFOAM skuznet OpenFOAM Running, Solving & CFD 99 March 16, 2017 06:07
problem with Min/max rho tH3f0rC3 OpenFOAM 7 February 23, 2013 06:37
Fluent msh and cyclic boundary cfdengineering OpenFOAM Other Meshers: ICEM, Star, Ansys, Pointwise, GridPro, Ansa, ... 48 January 25, 2013 04:28
Problem with rhoSimpleFoam matteo_gautero OpenFOAM Running, Solving & CFD 0 February 28, 2008 07:51
Trimmed cell and embedded refinement mesh conversion issues 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.