CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Gradient at the boundary

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 19, 2010, 06:09
Default Gradient at the boundary
  #1
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 208
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Hi all,

one important thing came to my mind yesterday and I would appreciate your comments on this! Let's say we have a fixedGradient BC, and consider evaluate method for our field:

Code:
template<class Type>
void fixedGradientFvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
    if (!this->updated())
    {
        this->updateCoeffs();
    }

    Field<Type>::operator=
    (
        this->patchInternalField() + gradient_/this->patch().deltaCoeffs()
    );

    fvPatchField<Type>::evaluate();
}
It means, that despite our selected scheme for the gradient calculation, we just use linear interpolation from the near-the-boundary cell center to the boundary face:

fi(wall) = fi(cell_center) + Grad0 * d

Is this a really correct estimation??? For example, if we use fourth order inside the domain for the interpolation and we end up with only (2nd?) order at the boundary using such BC type?

I would really appreciate your comments!

Have a nice day!

Alexander
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at

Last edited by makaveli_lcf; November 19, 2010 at 08:18.
makaveli_lcf is offline   Reply With Quote

Old   November 19, 2010, 09:01
Default
  #2
Senior Member
 
Francois
Join Date: Jun 2010
Location: Netherlands
Posts: 101
Rep Power: 6
Fransje is on a distinguished road
Dear Alexander,

The lower order estimation on the boundary is correct. Of course you would also like a higher-order estimation on the boundary, but because it makes the implementation of such a scheme more complex, a lower order estimation on the boundary is acceptable. It should not influence the order of convergence of the solution in your domain.
So your overall solution will not degrade from forth order accurate to second order accurate because you only use a second order accurate gradient estimation scheme on the boundary.

I hope this helps!

Kind regards,

Francois.
Fransje is offline   Reply With Quote

Old   November 19, 2010, 09:08
Default
  #3
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 208
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Francois, thank you for your reply!

Another point to discuss: there is no non-orthogonal correction at the boundary:

Code:
00161             //- Return gradient at boundary
00162             virtual tmp<Field<Type> > snGrad() const
00163             {
00164                 return gradient_;
00165             }
So, if we calculate boundary values using vector between the cell center and and the face center it is not correct especially for the sqewed cells.
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   November 21, 2010, 12:59
Default
  #4
Senior Member
 
Francois
Join Date: Jun 2010
Location: Netherlands
Posts: 101
Rep Power: 6
Fransje is on a distinguished road
Dear Alexander,

Could you not try a dot product between the velocity component and the the face area? Something like:
Code:
boundaryUxValue = ( Ux & yourBoundaryPatch.patch().Sf() );
etc.
Kind regards,

François.
Fransje is offline   Reply With Quote

Old   November 21, 2010, 15:47
Default
  #5
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 208
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
François, what should I check with this? Do you mean dot product between total velocity and normal vector to estimate flux?
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   November 22, 2010, 07:57
Default
  #6
Senior Member
 
Francois
Join Date: Jun 2010
Location: Netherlands
Posts: 101
Rep Power: 6
Fransje is on a distinguished road
Dear Alexander,

Yes. If you are afraid that the inbuilt implementations do not compute/use the solutions you expect, you can easily check that the values OpenFOAM uses are similar to those you would use in your own implementation. I thought this might help for your boundary.

Kind regards,

Francois.
Fransje is offline   Reply With Quote

Old   November 22, 2010, 11:20
Exclamation
  #7
Senior Member
 
Francois
Join Date: Jun 2010
Location: Netherlands
Posts: 101
Rep Power: 6
Fransje is on a distinguished road
Hum. A quick look at the wiki also returned the following answer. This should take care of the issues altogether.

Quote:
Here is a list of schemes using non-orthogonal correction on the boundary patches:
  • gradSchemes/fourthGrad
  • gradSchemes/extendedLeastSquaresGrad
  • gradSchemes/leastSquaresGrad
  • snGradSchemes/correctedSnGrad
  • surfaceInterpolation/schemes/skewCorrected
  • surfaceInterpolation/limitedSchemes/DeferredCorrectionLimitedScheme
  • surfaceInterpolation/limitedSchemes/LimitedScheme
  • surfaceInterpolation/limitedSchemes/linearUpwind

Kind regards,

Francois.
Fransje is offline   Reply With Quote

Old   December 27, 2010, 04:53
Default Correction at boundary
  #8
Member
 
Andrea Petronio
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 36
Rep Power: 7
andrea is on a distinguished road
Hi all,
can anyone clarify me this point, please:
chosen a corrected scheme for the Laplacian computation, the non-orthogonal correction is applied automatically to the boundaries too?

Thank you
Andrea
andrea is offline   Reply With Quote

Old   December 30, 2010, 11:05
Default Correction at boundary
  #9
Member
 
Andrea Petronio
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 36
Rep Power: 7
andrea is on a distinguished road
Hi all,
I'm confused about 2 points:

looking into http://foam.sourceforge.net/doc/Doxy...8C_source.html
I see on line 309:
Code:
// Boundary correction vectors set to zero for boundary patches
// and calculated consistently with internal corrections for
// coupled patches
so that the implementation of boundary correction described in Eugene PhD thesis
seems to be not implemented.

And I read from many source that DeltaCoeffs() should be simply the inverse
of the cell center to cell center vector while in the code is replaced by
Code:
// Stabilised form for bad meshes
DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
I don't know, maybe I'm looking in the wrong place!
Comments on these are appreciate
Andrea
andrea is offline   Reply With Quote

Old   January 13, 2012, 12:48
Default
  #10
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 535
Rep Power: 17
bigphil will become famous soon enough
Hi,

I have been looking at the calculation of gradients (fvc::grad) and surface-normal gradients (fvc::snGrad, and fvc::laplacian).

As far as I can see, the snGrad is calculated incorrectly on non-orthogonal fixedValue boundaries.
There should be non-orthogonal correction on fixedValue boundaries but it is set to zero (in surfaceInterpolation.C).

I have prepared a simple utility to check:
Code:
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fixedGradientFvPatchFields.H"

int main(int argc, char** argv)
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
  
  Info << "Testing field operations" << nl << endl;

  //-------------------------------------------------------------------------//
  //- define rotation
  //-------------------------------------------------------------------------//
  scalar theta = 1; // degrees
  const scalar PI = 3.141592;
  theta *= PI;
  theta /= 180;
  
  tensor rotMat(::cos(theta), -(::sin(theta)), 0,
		::sin(theta), ::cos(theta),    0,
		0,            0,               1);
  //-------------------------------------------------------------------------//



  //-------------------------------------------------------------------------//
  Info << "Reading alpha" << endl;

  volVectorField alpha
    (
     IOobject
     (
      "alpha",
      runTime.timeName(),
      mesh,
      IOobject::MUST_READ,
      IOobject::NO_WRITE
      ),
     mesh
     //     ((rotMat & mesh.C()) -  mesh.C()),
     //fixedValueCorrectedFvPatchVectorField::typeName
     //fixedValueFvPatchVectorField::typeName
     //fixedGradientCorrectedFvPatchVectorField::typeName
     //fixedGradientFvPatchVectorField::typeName
     //boundaryTypes
     );
  //-------------------------------------------------------------------------//



  //-------------------------------------------------------------------------//
  Info << nl << "Setting alpha field" << endl;
  //-------------------------------------------------------------------------//
  //- increment time
  runTime++;

  //- set the internal field
  alpha = ((rotMat & mesh.C()) -  mesh.C());

  //- set the fixedValue and fixedGradient boundaries
  forAll(alpha.boundaryField(), patchi)
    {
      if(alpha.boundaryField()[patchi].type() == fixedValueFvPatchVectorField::typeName)
	{
	  Info << "\tSetting alpha fixedValue boundaries for patch "
	       << mesh.boundary()[patchi].name() << endl;
	  const vectorField& patchC = mesh.boundaryMesh()[patchi].faceCentres();
      	  
	  alpha.boundaryField()[patchi] == ((rotMat & patchC) -  patchC);
      	}
      else if(alpha.boundaryField()[patchi].type() == fixedGradientFvPatchVectorField::typeName)
	{
	  Info << "\tSetting alpha fixedGradient boundaries for patch "
	       << mesh.boundary()[patchi].name() << endl;
	  vectorField n = mesh.boundary()[patchi].Sf() / mesh.boundary()[patchi].magSf();
      	  fixedGradientFvPatchVectorField& alphaPatchGrad = refCast<fixedGradientFvPatchVectorField>(alpha.boundaryField()[patchi]);
      	  
	  //- analytical internal gradAlpha is (rotMat.T() - I)
	  //- therefore analytical boundary gradient is:
	  alphaPatchGrad.gradient() = n & (rotMat.T() - I);
      	}
      else
	{
	  FatalError << "The geometry should only have boundaries of type"
		     << " fixedValue or fixedGradient!"
		     << exit(FatalError);
	}
    }

  alpha.write();
  //-------------------------------------------------------------------------//


  surfaceVectorField n = mesh.Sf()/mesh.magSf();


  //-------------------------------------------------------------------------//
  //- gradient of vector field
  //- least squares seems right in the boundary cells...
  //- Gauss linear is wrong in all skewed cells as expected
  //-------------------------------------------------------------------------//
  Info << nl << "Writing out the gradient" << endl;

  volTensorField gradAlpha
    (
     IOobject
     (
      "gradAlpha",
      runTime.timeName(),
      mesh,
      IOobject::NO_READ,
      IOobject::NO_WRITE
      ),
     fvc::grad(alpha)
     );
  Info << "\tMax mag is " << max(mag(gradAlpha)).value()
       << "\tMin mag is " << min(mag(gradAlpha)).value()
       << endl;
    
  gradAlpha.write();
 //-------------------------------------------------------------------------//
 

  //-------------------------------------------------------------------------//
  //- laplacian of alpha
  //- lapAlpha is fine when the BC of alpha are fixedValue but lapAlpha is
  //- wrong in the boundary cells when the BC is fixedGradient... assuming I
  //- set the fixedGradient BCs correctly
  //-------------------------------------------------------------------------//
  Info << nl << "Writing out the laplacian of alpha" << endl;
  volVectorField lapAlpha
    (
     IOobject
     (
      "lap(alpha)",
      runTime.timeName(),
      mesh,
      IOobject::NO_READ,
      IOobject::NO_WRITE
      ),
     fvc::laplacian(alpha, "lap(alpha)")
     // equivalent to fvc::surfaceIntegrate(fvc::snGrad(alpha)*mesh.magSf())
     );
  Info << "\tMax mag is " << max(mag(lapAlpha)).value()
       << "\tMin mag is " << min(mag(lapAlpha)).value()
       << endl;
  lapAlpha.write();
  //-------------------------------------------------------------------------//


  Info << nl << endl;
  return 0;
}
This utility defines an analytical vector field alpha.
grad(alpha) should have the constant components (-0.00015, 0.017452, 0, -0.017452, -0.00015, 0, 0, 0, 0) and laplacian(alpha) should be (0, 0, 0).The alpha field actually corresponds to rigid rotation about the z-axis.

fvc::grad(alpha) is calculated and seems correct, but the fvc::laplacian(alpha) is incorrect on non-orthogonal fixedValue boundaries!

I have attached a non-orthogonal test case which shows the problem.

Does anybody have any thoughts on this, am I misunderstanding something?

Philip


Edit: by the way, fvc::laplacian(alpha) is correct for fixedGradient boundaries because they don't need non-orthogonal correction as the surface-normal gradient is already known.

Edit2: I tried this on OpenFOAM-2.1.x and OpenFOAM-1.6-ext with the same results. Also I checked old foam code and it seems to have no correction on the boundaries too, hmnnn...
Attached Images
File Type: jpg Laplacian incorrect on boundaries.jpg (75.2 KB, 25 views)
Attached Files
File Type: zip nonOrthFixedValue.zip (82.5 KB, 4 views)
bigphil is offline   Reply With Quote

Old   November 25, 2012, 15:26
Default
  #11
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 535
Rep Power: 17
bigphil will become famous soon enough
Hi,

for future reference,

OpenFOAM sets non-orthogonal correction to zero on the boundaries. As described in Hrv's thesis, it is assumed that the dependent variable is uniform along the boundary face and hence no non-orthogonal correction is applied.

However, there are cases (for instance in solid mechanics) when this is a poor assumption and it is assumed that the dependent variable varies linearly along the boundary face and boundary non-orthogonal corrections are imperative.

So if you need boundary non-orthogonal correction then custom boundary conditions must be used where non-orthogonal corrections are imposed (see here).

Philip
bigphil 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
External Radiation Boundary Condition for Grid Interface CFD XUE FLUENT 0 July 9, 2010 02:53
External Radiation Boundary Condition (Two sided wall), Grid Interface CFD XUE FLUENT 0 July 8, 2010 06:49
mesh file for flow over a circular cylinder Ardalan Main CFD Forum 6 April 17, 2010 23:40
normal temperature gradient on a boundary Sandrine FLUENT 0 June 10, 2009 12:29
Boundary condition vector field gradient from two sides of patch face quba OpenFOAM 0 December 12, 2007 05:26


All times are GMT -4. The time now is 19:55.