[alberto]

After thinking about it a bit more, I tend to think what is done in OpenFOAM is correct, even if probably not the expected behaviour. I see a boundary condition as an information specified by the user, and the code should not alter it, if the user does not allow that.

I think what is confusing here is that in the "textbook" approach, we actually mix the discretization of the equation and of the boundary to implicitly include this effect. In OpenFOAM the two aspects are separate: the div scheme does its job of discretizing the terms of the equation, while the BC takes care of the value at the boundary, which, for a fixedValue, has nothing to do with how you perform the interpolation, if you think to this condition as something that strictly applies a value on the boundary faces.

In your case, as Hrv said, you need a condition where the value specified on the face is not enforced explicitly, as fixedValue does, but it is used to define a gradient, based on the internal properties of the flow, so that T on the face becomes the specified value. You'll have to solve this implicitly.

In this way, however, there is no need to change a significant portion of the code, and it won't affect other boundary conditions.

[eugene]

There is one fundamental problem with your proposal, it does not distinguish between the gaussConvection and gaussLaplacian operators. The really fundamental aspect of what Santiago is proposing is that the boundary can have a different value depending on the interpolation scheme of the specific operator. Thus, for a fixed value outlet, the outlet face for transported property T would have the following values depending on the operator:

div(phi,T) Gauss upwind; T_f = T_i
laplacian(kappaEff,T) Gauss linear uncorrected; T_f = T_b

where T_i is the internal value, T_b is the boundary value and T_f is the value used to construct the matrix. (For the laplacian it would actually be a gradient not a value, but the principle is the same.) It is not possible with the current way in which weights are constructed to distinguish between upwind and linear interpolation at the face. To achieve the distinction requires a straight-forward, but non-trivial alteration of all valueInternalCoeffs and valueBoundaryCoeffs functions and also the weights function in fvPatch.C

This does not mean that the basic boundary condition behaviour has to be altered. It would still be possible to maintain equivalence to the current scheme and add new boundaries like: fixedValueSantiago and zeroGradientSantiago that behave in a way consistent with Santiago's proposal.

[alberto]

Let's assume we implement your approach. What value do you assign to the face when you store the solution? You end up with two different values at the same face...

Maybe I am a bit stubborn :-), but that is why I think the currently available approach is fine. It should be the boundary to provide a consistent value (unique) to the operators, more than the operator to use different face values.

[santiagomarquezd]

Eugene's comment is a correct summation of my proposal (really is not my proposal but what I learnt from my advisor). The value written in the solution is the same that user gave us, because it is fixed (the tricks are done only internally). The actual implementation is consistent in the idea of face value given by user is the boss in all the schemes, but this concept leads to unstable numerical schemes in the case of upwind. I'm a high enthusiast of FOAM but in this case we have to assume that even when actual machinery is very robust this concept is wrong, or at least I haven't seen a better explanation.
Until now, as Eugene commented, the problem is that is necessary to pass some other data to BC class in order it can manage different div and grad schemes. We've tried to solve the problem too but we finally figured that Eugene said. Another path is to add some attributes to BC class and initialize them in the BC constructor, for example, reading the scheme from fvSchemes, then when internalCoeffs are calculated, upwind particular case can be taken in account. The problem is that many fvm::div operators can be assembled in the same set of equations and it is not possible to know on the fly which of these terms is being assembled.

Here we keep trying.

Regards.

[eugene]

Santiago, well it seems to me that a simple switch would be half a job. If your div scheme uses 50/50 upwind+linear blending, then your face value should consist of 50% internal value + 50% face value, no? The machinery for passing the required information to valueInternalCoeffs and valueBoundaryCoeffs is already in place, it is just not being used properly.

Currently, all differencing schemes assign a weighting to the boundary value, irrespective of what the interpolation scheme is. In no instance other than the coupled boundary is this weighting actually used. For most non-coupled boundaries, the weighting is simply set to 1 and ignored. We can leverage this weighting interface designed to set the interpolation on coupled boundaries to set the interpolation on non-coupled boundaries as well. Since the actual value of the weighting is not used by the existing non-coupled boundaries, messing around with its values will not affect any of the current implementations.

[alberto]

Call me stubborn , but I am not convinced of the need of messing with the interpolation schemes to solve this problem. I'll try to explain.

Let's examine first what happens physically in your system (even if your problem is indeed a bit far from reality), since it helps to set conditions properly. You have a flow with constant temperature T0, and you convect it against a boundary at a constant, but different temperature T1. To have this, you assume a localized, instantaneous change in temperature (not likely to happen in reality).
However, mathematically this is a possible problem, and the numerical instability that arises from the current implementation is that the boundary condition is not appropriate (the scheme has nothing to do with it, it is just a case that with some non-biased scheme it works).

Your boundary condition strongly enforces a value at the face, independently from what is convected against it, while a proper implementation would impose the boundary value more weakly (it is really a common strategy in this type of problems), by ensuring that the proper heat exchange takes place. This corresponds to enforcing the gradient, and not the value, which will result as a consequence.

If you do that, you respect the physics of the problem, and you do not have to play any trick with the numerics at all, avoiding inconsistencies in the values used in the discretization of your equation.

Best,

[eugene]

Here is my suggestion for a simple weighting scheme on the boundary:
pure upwind = 1
pure downwind = 0
central/linear = 0.5

Since we don't know anything concrete about what is outside the boundary, lets limit the weighting such that the boundary will not pretend to know anything more than what the user provided value and the internal domain can provide. This means that for inlet flows, the weighting can never be larger than 0.5 and for outlet flows, the weighting can never be smaller than 0.5.

The actual face weighting value before limiting should just be the relative contribution of the owner cell to the face value. We have assumed that using the specified face value corresponds to central weighting of 0.5 (0.5 owner, 0.5 neighbour). If we use a blend of 50% upwind + 50% linear, this equates to a owner weighting of 0.75. So to convert our weighting factor into an equation to combine the internal value and the face value for fixed value outlets:

valueInternalCoeffs = max(0, (w - 0.5)/0.5)

valueBoundaryCoeffs = min(1, (1 - w)/0.5) * (*this)

The limiting ensures that we never use more downwind than upwind interpolation at the outlet and generally apply the fixed boundary value at the inlet (unless you specify a downwind scheme!). Note, that under no circumstances will the boundary "know" more about what happens outside the domain than supplied by the fixed value boundary.

[santiagomarquezd]

Originally Posted by alberto

To have this, you assume a localized, instantaneous change in temperature (not likely to happen in reality).

I can't understand this point, solution has a big gradient at the outlet, but this gradient is, in fact, which allows to remove the magnitude by means of the difussion at outlet. At the inlet the flux is almost purely advective, at the outlet is the opposite. This problem appears due the kind of discretization obtained by cell centered FVM. Using FEM or FDM you can use a node in the boundary so the first and last equations are simply T_1=T_BC1, T_n=T_BC2, and the upwind is applied to all other nodes without problems. Using this discretizations the solution is practically T_(n-1)~T_BC1 and T_n=T_BC2. There is a jump too, but it is due the discretization again.

Originally Posted by alberto

This corresponds to enforcing the gradient, and not the value, which will result as a consequence.

Could you put this idea in mathematical notation, really I desire to understand, but I can't do it.

Regards.

[eugene]

The thing is Alberto, there are a wide range of (admittedly marginal) problems that will be better served by, lets call it a scheme dependent boundary. In every case, you could go and write a specific boundary that might alleviate the problem as you propose, but why not simply provide this option (that incidentally should reproduce text-book results). Also, I have several experiences where specifying a gradient as you propose does not guarantee stability. Granted, most people will never encounter these situation (convection of vector fields other than velocity), but it doesn't change the fact that there are many conditions where a scheme dependent boundary can be advantageous because of its generality.

In any case, the actual amount of "messing" with the code is quite modest and as far as I can tell and can be done in such a way that it has no effect on the existing boundary schema.

Let me give you another example, a colleague of mine was doing some moving ground aerodynamic simulations. Because of the vagaries of snappy, some of the ground boundary faces near the wheels were not perpendicular to the direction of motion leading to a net flux into some of the ground boundary faces. Since the these "outlets" are "fixed value", they cause all kinds of problems, even with a near-linear interpolation scheme. Yes, we have more advanced boundaries that were built specifically to deal with this kind of situation, but the scheme dependent boundary would be able to handle this in its stride and many similar problems besides.

[santiagomarquezd]

Originally Posted by eugene

Currently, all differencing schemes assign a weighting to the boundary value, irrespective of what the interpolation scheme is. In no instance other than the coupled boundary is this weighting actually used. For most non-coupled boundaries, the weighting is simply set to 1 and ignored.

I was checking this today and for upwind the weights at boundaries are set in the same form that internal ones. 1 for downwind and 0 for upwind (they are set in 1 in the constructor but are changed afterwards). So an easy way to hack the fvm::div method is doing valueInternalCoeffs=weights if scheme=upwind and BC=fixedValue. One implementation is to do it directly in gaussConvectionScheme.C by an if statement, another is, I think, the way of Eugene in the last post, but I have to study it to understand it at all.

Regards.

[alberto]

Originally Posted by santiagomarquezd

I can't understand this point, solution has a big gradient at the outlet, but this gradient is, in fact, which allows to remove the magnitude by means of the difussion at outlet.

Yes, physically this is what happens. But the numerical setup that gives you problems does not say exactly that, due to the upwind interpolation.

Could you put this idea in mathematical notation, really I desire to understand, but I can't do it.

Just think to a condition where the imposed gradient is dynamically determined from the upstream value and the boundary value.

[santiagomarquezd]

Eugene,

Originally Posted by eugene

Here is my suggestion for a simple weighting scheme on the boundary:
pure upwind = 1
pure downwind = 0
central/linear = 0.5

Since we don't know anything concrete about what is outside the boundary, lets limit the weighting such that the boundary will not pretend to know anything more than what the user provided value and the internal domain can provide. This means that for inlet flows, the weighting can never be larger than 0.5 and for outlet flows, the weighting can never be smaller than 0.5.

The actual face weighting value before limiting should just be the relative contribution of the owner cell to the face value. We have assumed that using the specified face value corresponds to central weighting of 0.5 (0.5 owner, 0.5 neighbour). If we use a blend of 50% upwind + 50% linear, this equates to a owner weighting of 0.75. So to convert our weighting factor into an equation to combine the internal value and the face value for fixed value outlets:

valueInternalCoeffs = max(0, (w - 0.5)/0.5)

valueBoundaryCoeffs = min(1, (1 - w)/0.5) * (*this)

The limiting ensures that we never use more downwind than upwind interpolation at the outlet and generally apply the fixed boundary value at the inlet (unless you specify a downwind scheme!). Note, that under no circumstances will the boundary "know" more about what happens outside the domain than supplied by the fixed value boundary.

if I understand well your proposal is to use the weights passed to valueInternalCoeffs and valueBoundaryCoeffs as a way to know what scheme is being used, so 1 is used for a face downstream of the last cell center, 0 for the opposite, both in upwind scheme, and 0.5 for centered/linear scheme. This would be a slight difference respect of actual weights due central/linear are 1 in stead of 0.5. Is my intepretation correct? If it were this way, how about other schemes like QUICK, TVD, NVD? I think they are performing well now, so it would be enough with giving them values not equal to those are used for upwind. Is that correct?

Regards.

[eugene]

Yes. The formulation I suggested will result in the correct proportional combination of upwind and the boundary value for any scheme that is a combination of upwind and downwind provided the input owner weighting was correct. To get the correct behaviour, you would however have to modify pretty much every non-upwind interpolation scheme you want to use such that it provides the correct owner weighting on non-coupled boundaries. Your main problem currently is that the coupled boundaries rely on the weighting for central schemes being 1, so you cannot simply reuse the current coupled boundary weighting functions without alteration on non-coupled boundaries.

But you don't need to worry about all that to begin with. To start, you just need to consider the upwind and linear schemes. For upwind the boundary weightings are already correct. For linear, you need to edit the weights function in linear.H such that it returns 0.5 on the non-coupled boundaries instead of the current 1. Then all you need to do is create a derivative of fixed value boundary that formulates the valueInternal and valueBoundaryCoeffs in the way I described and you should be done.

[Duncan_vdH]

I am following this discussion with interest but I do not get the idea.
We actually use the specific case of 1-D steady convection-diffusion equation with two-sided dirichlet bc's to show that these bc's lead to a well posed problem in a course I am teaching.

However, why would you want to interpolate the scalar variable, if the value you need for the flux is present where you need it?
If you use this approach, so use this exact value for the flux (no interpolation) you get the correct numerical solution with first order upwind (little bit overestimated boundary layer thickness) or with a central scheme (when the mesh peclet number < 2) on a cell-centered grid.

So I am really surprised that OF would give the wrong solution.

[santiagomarquezd]

Eugene, I think your solution is the most elegant, I'll try to implement it. BTW I'll look for some FVM examples of this problem, because all my references are in FEM and FDM.

Duncan, I attached a FOAM case at the beggining of the thread, check it if you wish. Problem appears with upwind in outlet BC due the using of downstream value in a upwinded scheme.

Regards.

[Duncan_vdH]

Hmmm, I seem to get the wrong solution even with a central discretisation.
What settings did you use to get the correct answer (with the correct boundary layer)??
The steadyState seems to be ignored...

Regards,
Duncan

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;
div(phi,T) Gauss limitedLinear 1;
}

laplacianSchemes
{
default none;
laplacian(DT,T) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
T ;
}

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 0 0 0 0 0 0];

internalField uniform 2;

boundaryField
{
frontAndBackPlanes
{
type empty;
}
Pared
{
type zeroGradient;
}
Entrada
{
type fixedValue;
value uniform 1;
}
Salida
{
type fixedValue;
value uniform 2;
}
}

[santiagomarquezd]

Hi, the way to assemble fvm::div operator that we were discussing is explained in Versteeg (1995) book p. 116. There he uses the flux at the boundary, but the value at upwind cell centroid not the boundary value. Thanks to Emanuele for giving the reference (http://www.cfd-online.com/Forums/ope...tml#post280021).

Regards

[eugene]

Well that seals it in my opinion.
Good timing too, I was just wondering about how you were getting along as I have run into another problem where this div thing at the outlet is an issue. I will hack around a bit and let you know if I get any interesting results.

[santiagomarquezd]

Eugene, thx for answer. My advisor's FOAM Matlab emulator is now kicking and running, at least for the scalarTransportFoam emulation. We're running some cases in which we've found more problems with advection when they are run in FOAM with upwind. Once my advisor finished the runnings I'll run the same cases but with the patched version of FOAM (this version is not so elegant as as your solution but does the job), if there would be some issues related this thread I'll post them.

Let me know of your findings too.

Regards.

[eugene]

The correct formulation turned out to be unnecessary for me as there was no flux over the new boundaries I had to introduce.