[hjasak]

Why don't you just write a direction mixed boundary condition which has got a separate refValue, refGrad and valueFraction for the convection and diffusion part? Each operator does its own matrix assembly anyway and you have 4 separate functions for the coefficient:

virtual tmp<Field<Type> > valueInternalCoeffs
virtual tmp<Field<Type> > valueBoundaryCoeffs

virtual tmp<Field<Type> > gradientInternalCoeffs() const
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const

Hrv

[jorkolino]

Thanks Santiago for the help. I have a debugger, I only need to figure out where to supply the above commands.
Thanks Hrv for your help too, I have to figure out for myself the meaning of this .

Could anybody please provide little more info on how to use the lines suggested by hrv? These are virtual templated functions? Lets say I have a patch called "outlet" and a variable called LMA. What should I do to enquire the coefficients (and possibly change them) in the near outlet cells of the variable LMA? Anybody ?

[jorkolino]

Originally Posted by alberto

If you set DT = 0.001, with that case, without relaxing, refining the grid next to the outlet and/or using a small time stepping in the case of unsteady runs, you won't get the correct solution for the reason explained above (you lose diagonal dominance).

The pictures attached show you what I get with OF 1.7.x, using scalarTransportFoam with

T0 = 1
T1 = 2
DT = 0.001
Upwind

and with the same settings, but relaxing with URF 0.001 (still wrong solution), and 0.00001. The case and code are attached, since I had to modify scalarTransportFoam to include relaxation.

P.S. Plots are made in paraview, so the scales are what they are...

Hi there,

What is the correct precedence in the code (cause I've seen it both ways): first solve then relax or vice versa?
I noticed relaxation factor of unity has huge impact on results when linear div scheme is used, in comparison to when no relaxation is used. I would expect relaxation factor of 1 to have no impact on the solution. That happens to a well defined parameter: fixed inlet value and zerogradient outlet value in 1D.
It is almost certain that there is a problem with Foam not only with certain BC, but also with the precedence of relaxation. URF of 1 should not change the behavior at all.

[dhuckaby]

George,

from your discussion, its seems for your equation you are trying to advect information from the outlet into the domain using the negative of the velocity. As such it would seem that you would want to use "-phi" within your divergence operator, e.g. -div(-phi, LMRA).

Dave

[alberto]

Originally Posted by jorkolino

Hi there,

What is the correct precedence in the code (cause I've seen it both ways): first solve then relax or vice versa?

Relaxation is applied implicitly before solving if you do:

yourEqn = ...;
yourEqn.relax();
yourEqn.solve();

If you see the relax() method applied to a field (like p in simple algorithm), you are explicitly relaxing the variable.

You typically relax the equation before solving the equation. In the SIMPLE algorithm, you explicitly relax the pressure after the solution of the pressure equation and the flux correction, but before the velocity correction.

I noticed relaxation factor of unity has huge impact on results when linear div scheme is used, in comparison to when no relaxation is used. I would expect relaxation factor of 1 to have no impact on the solution. That happens to a well defined parameter: fixed inlet value and zerogradient outlet value in 1D.
It is almost certain that there is a problem with Foam not only with certain BC, but also with the precedence of relaxation. URF of 1 should not change the behavior at all.

There is no problem in the relaxation itself, which, if done properly, should not change the converged solution at all in general (if you exclude URF = 0). There is some sensitivity to under-relaxation in OpenFOAM due to the co-located grid arrangement, however, as discussed in another thread it is usually smaller than the discretization error of the corresponding schemes (see also Peric's book for the theory).
What you see is due to something else (can your solution converge with URF = 1?)

About the directioMixed BC, it is basically an extended Robin condition (a*U + b*dU = c), which allows you to specify a different "fraction" of fixedGradient / fixedValue in each direction (Sorry for the horrible wording ).

There are a couple of BC's to look at, which are derived from the directionMixed BC:

- outletStabilised, which applies upwind differencing to outlets
- pressureInletOutletVelocity (see header for description)

Best,

[jorkolino]

Originally Posted by hjasak

Why don't you just write a direction mixed boundary condition which has got a separate refValue, refGrad and valueFraction for the convection and diffusion part? Each operator does its own matrix assembly anyway and you have 4 separate functions for the coefficient:

virtual tmp<Field<Type> > valueInternalCoeffs
virtual tmp<Field<Type> > valueBoundaryCoeffs

virtual tmp<Field<Type> > gradientInternalCoeffs() const
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const

Hrv

Hi,

I appreciate your help. If you could please briefly outline the steps I need to follow. I imagine it should be like this:

1. First I define a new class inheriting from fvFixedValuePatch, similar to the example oscillatingFixedValueFvPatchField.H (http://foam.sourceforge.net/doc/Doxy...8H_source.html).
2. Then I need to define my inlets(outlets) in the dictionary file to reflect this new BC.
3. I should print(and modify later on) the 4 contributions to the coefficients using the 4 functions above.

Is it correct? Thanks for your help.

[hjasak]

Nope, but you are on the right track. Shall we make a deal - I will help you do this and we'll then include it into -Extend?

This is what you do:
- take a mixed condition and make a new one, call it doubleMixed or valueGradMixed. Be creative with the name, I don't mind
- we need two copies of:

//- Value field
Field<Type> refValue_;

//- Normal gradient field
Field<Type> refGrad_;

//- Fraction (0-1) of value used for boundary condition
scalarField valueFraction_;


You can call them

valueRefValue valueRefGrad valueValueFraction
gradientRefValue gradientRefGrad gradientValueFraction

- pull them through constructors, mapping, access functions, autoMap and rMap,

- we have 2 functions which become a BIG problem; we will leave them for later. They are:

snGrad() and evaluate(). For the moment, use the value entries and we'll figure them out later.

Now you have 4 "impact" functions:


//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueInternalCoeffs
(
const tmp<scalarField>&
) const;

//- Return the matrix source coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueBoundaryCoeffs
(
const tmp<scalarField>&
) const;


AND

//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientInternalCoeffs() const;

//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const;

In the valueInternalCoeffs and valueBoundaryCoeffs you will use the valueRefValue valueRefGrad valueValueFraction and in gradientInternalCoeffs() and gradientBoundaryCoeffs() you will use the
gradientRefValue gradientRefGrad gradientValueFraction.

The formulae are the same, it's just that in the value and gradient entries you allow a separate mix of fixedValue to fixedGradient, because you have 6 controls now.

I think you've now got the drift; it remains to sort out the snGrad and evaluate, but I guess the genii who claim this is physical behaviour will tell me what to do.

Hrv

[jorkolino]

Originally Posted by alberto

What you see is due to something else (can your solution converge with URF = 1?)

It converges, slowly, but I am surprised as to why URF of 1 should change at all the matrix coefficients...I refer to Patankar on that matte, but obvioulsly I miss something here.

I calculate a S-S problem on a transient solver (scalarTransportFoam) with ddt set to SteadyState. It converges on several time steps rather slowly with URF=1, because it does a single iteration per time step on the calculated variable. Instead I would expect it to iterate until the final solution is reached on the first time step. Maybe I have to use explicit relaxation here.

Originally Posted by alberto

There are a couple of BC's to look at, which are derived from the directionMixed BC:

- outletStabilised, which applies upwind differencing to outlets
- pressureInletOutletVelocity (see header for description)

Best,

Thanks, I think I can make use of the first one. But I still need a method to remove the diffusive part at inlets for my second parameter, so I need to search for way to do that.I have asked two posts above wether my suggestion is correct.

[jorkolino]

Originally Posted by hjasak

Nope, but you are on the right track. Shall we make a deal - I will help you do this and we'll then include it into -Extend?

Whatever this means, I will be proud to

Originally Posted by hjasak

I think you've now got the drift; it remains to sort out the snGrad and evaluate, but I guess the genii who claim this is physical behaviour will tell me what to do.

Hrv

Thanks for the credit but the roots of this parameter go back to Ad Roos and his thesis. It is a mere mathematical equation with BC available for all boundaries and as such it should be possible to obtain unique solution. Question is, can we do it with control volume method? I should think it can be done.

l attach the header and source files, together with Make files to create a libraby. Now everything compiles well. I've used dummy expressions in the snGrad and evaluate functions, just to pass the compile stage. I still need your help on how to split the effect of value and grad in those functions.

[jorkolino]

Originally Posted by dhuckaby

George,

from your discussion, its seems for your equation you are trying to advect information from the outlet into the domain using the negative of the velocity. As such it would seem that you would want to use "-phi" within your divergence operator, e.g. -div(-phi, LMRA).

Dave

Hi Dave,

thanks for reply. Mathematically the expression above seems to me equivalent to div(phi,LMRA), which is originally used. But I may try your suggestion anyway. The big problem now is how to set up B.C. such that to ensure stable matrix, and also how to disable diffusive contribution at inlets when using fixedValue B.C at inlets.

Regards

[jorkolino]

Hi there,

I am posting the boundary condition which uses upstream value for the representation of the convection term. I have tested it for upwind scheme and can confirm that solution stays bounded for any value of Peclet. On the contrary, the basic fixedValue B.C. oscillated under same conditions when Pe>2.

The problem that I solve requires that the outlet value is fixed, but the inlet is unknown. When using zeroGradient for inlet I have problems with both the standard fixedValue and the atttached B.C. They create instability in case when source terms are present. So, I was thinking on the following.

After integration of the PDE for the last cell, there is no mathematical restriction on how many B.C will be defined for that last cell. In other words, I can apply fixed value and fixed gradient on the resulting integral equation, thus creating two discrete equations. So, why not generate two equations for the near-boundary cell, and skip an equation for the, lets say, the first cell? The matrix will be solvable, and the values in all cells will emerge as far as No of equations is equal to No of cells. Any ideas how can this be done with OF?