CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Bugs

Wrong fvm::div assembling

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

Like Tree17Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 20, 2010, 08:12
Default
  #81
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   December 20, 2010, 09:29
Default
  #82
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
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 ?

Last edited by jorkolino; December 21, 2010 at 11:40.
jorkolino is offline   Reply With Quote

Old   December 20, 2010, 13:44
Default
  #83
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
Quote:
Originally Posted by alberto View Post
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.

Last edited by jorkolino; December 21, 2010 at 03:53.
jorkolino is offline   Reply With Quote

Old   December 21, 2010, 14:24
Default
  #84
New Member
 
David Huckaby
Join Date: Jul 2009
Posts: 21
Rep Power: 7
dhuckaby is on a distinguished road
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
dhuckaby is offline   Reply With Quote

Old   December 21, 2010, 15:58
Default
  #85
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by jorkolino View Post
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.

Quote:
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,
fumiya likes this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; December 21, 2010 at 15:59. Reason: clarified
alberto is offline   Reply With Quote

Old   December 22, 2010, 05:46
Default
  #86
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
Quote:
Originally Posted by hjasak View Post
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.
jorkolino is offline   Reply With Quote

Old   December 22, 2010, 05:59
Default
  #87
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   December 22, 2010, 06:19
Default
  #88
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
Quote:
Originally Posted by alberto View Post
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.

Quote:
Originally Posted by alberto View Post
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 is offline   Reply With Quote

Old   December 22, 2010, 10:08
Default
  #89
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
Quote:
Originally Posted by hjasak View Post
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


Quote:
Originally Posted by hjasak View Post

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.
Attached Files
File Type: gz mixedValueGrad.tar.gz (3.2 KB, 14 views)

Last edited by jorkolino; December 28, 2010 at 04:26.
jorkolino is offline   Reply With Quote

Old   December 22, 2010, 12:11
Default
  #90
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
Quote:
Originally Posted by dhuckaby View Post
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 is offline   Reply With Quote

Old   December 27, 2010, 13:54
Default
  #91
Member
 
George Pichurov
Join Date: Jul 2010
Posts: 39
Rep Power: 6
jorkolino is on a distinguished road
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?
Attached Files
File Type: gz fixedValueZeroConvFlux.tar.gz (5.3 KB, 21 views)

Last edited by jorkolino; January 19, 2011 at 04:50.
jorkolino is offline   Reply With Quote

Reply

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
Warning: Dynamic zone with wrong CG using 6DOF Manoj Kumar FLUENT 1 August 11, 2012 04:03
BuoyantBoussinesqSimpleFoam and axial-symmetric 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


All times are GMT -4. The time now is 18:34.