
[Sponsors] 
September 27, 2010, 15:04 

#21 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
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. 

September 28, 2010, 08:55 

#22 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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 straightforward, but nontrivial 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. Last edited by eugene; September 28, 2010 at 09:07. Reason: added another paragraph 

September 29, 2010, 09:37 

#23 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
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. 

September 29, 2010, 10:33 

#24 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar Last edited by santiagomarquezd; September 29, 2010 at 10:36. Reason: Spelling mistakes 

September 29, 2010, 12:27 

#25 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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 noncoupled 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 noncoupled boundaries as well. Since the actual value of the weighting is not used by the existing noncoupled boundaries, messing around with its values will not affect any of the current implementations. 

September 29, 2010, 14:30 

#26 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
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 nonbiased 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, 

September 29, 2010, 17:20 

#27 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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. 

September 29, 2010, 17:29 

#28  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
Quote:
Quote:
Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

September 29, 2010, 17:43 

#29 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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 textbook 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 nearlinear 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. 

September 29, 2010, 17:45 

#30  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
Quote:
Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

September 29, 2010, 18:30 

#31  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
Quote:


September 30, 2010, 19:57 

#32  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
Eugene,
Quote:
Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

October 1, 2010, 10:50 

#33 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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 nonupwind interpolation scheme you want to use such that it provides the correct owner weighting on noncoupled 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 noncoupled 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 noncoupled 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. 

October 4, 2010, 15:38 

#34 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7 
I am following this discussion with interest but I do not get the idea.
We actually use the specific case of 1D steady convectiondiffusion equation with twosided 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 cellcentered grid. So I am really surprised that OF would give the wrong solution. 

October 4, 2010, 16:29 

#35 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

October 5, 2010, 10:42 

#36 
New Member
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7 
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; } } 

October 20, 2010, 13:47 

#37 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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 (Convectiondiffusion in 1D : wrong solution for a large Delta x).
Regards
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

October 21, 2010, 04:55 

#38 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
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. 

October 22, 2010, 17:43 

#39 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

October 25, 2010, 03:42 

#40 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12 
The correct formulation turned out to be unnecessary for me as there was no flux over the new boundaries I had to introduce.


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Warning: Dynamic zone with wrong CG using 6DOF  Manoj Kumar  FLUENT  1  August 11, 2012 04:03 
BuoyantBoussinesqSimpleFoam and axialsymmetric results wrong mass flow  Thomas Baumann  OpenFOAM  6  December 21, 2009 11:31 
udf error  srihari  FLUENT  0  February 9, 2009 10:00 
Pressure contour seems wrong???  Harry Qiu  FLUENT  1  June 29, 2001 05:53 
What's wrong with my UDF?  olivia  FLUENT  1  June 23, 2001 17:06 