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

flux seems not conserved in my modified scalarTransportFoam

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

Reply
 
LinkBack Thread Tools Display Modes
Old   October 3, 2009, 08:59
Default flux seems not conserved in my modified scalarTransportFoam
  #1
New Member
 
Daniel Riemann
Join Date: Oct 2009
Posts: 3
Rep Power: 7
danielr is on a distinguished road
I try to understand why the flux conservation in a simple diffusion problem seems to be inaccurate with my following implementation:

I have modified the scalarTransportFoam solver to accept a nonuniform diffusivity k instead of the constant diffusivity DT. The velocity field U is set zero. So I solve the diffusion problem
d/dt T=div k grad T
in a unit square (16x16 cells) with k=0.01 in the upper half and k=1 in the lower half.
Boundary conditions are T=0 at the top and fixedGradient at the bottom.

The attached figure shows a color plot of T (above) at steady state (integration time t=100) and a probe line (below) of k,T and the flux magnitude |F| along the vertical axis.

I calculate the diffusive flux (F=k grad T) in the modified scalarTransportFoam as F=k*fvc::grad(T).
The calculated F shows a wiggle, where k jumps from 1 to 0.01, which can't be correct in steady state.

What is the correct expression to evaluate the flux k*grad T?

Thanks
Daniel
Attached Images
File Type: jpg flux.jpg (35.6 KB, 48 views)
danielr is offline   Reply With Quote

Old   October 3, 2009, 17:50
Default
  #2
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
Hi Daniel,

When you solve a laplacian with a finite volume method, the fluxes you use are evaluated on cell *faces*. The flux you computed from a gradient in cell centres (using an fvc operator) therefore does'n have anything to do with the fluxes you actually use in the solution. So, I'm not surprised it looks strange. The fluxes across cell faces are certainly conserved in the code.

Another thing: When you have such great variation in the diffusivity, you should be careful with the interpolation scheme you choose for k. It's physically more sensible to use "harmonic" rather than "linear" interpolation (especially if you have interfaces where k has step changes). Remember that a laplacian is evaluated as a sum of all face fluxes k_f*snGrad(T), where k_f is the interpolated value of k at each face.

When you select a laplacian scheme like "Gauss linear corrected" in fvSchemes, then "Gauss" stands for this method of surface integrals on faces (only option, I think), "linear" is the interpolation scheme for k in your case (here I would use "harmonic"), and "corrected" specifies that the surface-normal gradient of T is computed with full non-orthogonal corrections applied.

/Ola
olwi is offline   Reply With Quote

Old   October 4, 2009, 12:39
Default
  #3
New Member
 
Daniel Riemann
Join Date: Oct 2009
Posts: 3
Rep Power: 7
danielr is on a distinguished road
Quote:
Originally Posted by olwi View Post
Hi Daniel,

When you solve a laplacian with a finite volume method, the fluxes you use are evaluated on cell *faces*. The flux you computed from a gradient in cell centres (using an fvc operator) therefore does'n have anything to do with the fluxes you actually use in the solution. So, I'm not surprised it looks strange. The fluxes across cell faces are certainly conserved in the code.
I was suspecting something like this. I didn't doubt that the scalarTransportFoam was correct, but was not sure how to get the correct flux in my case.

Quote:
Remember that a laplacian is evaluated as a sum of all face fluxes k_f*snGrad(T), where k_f is the interpolated value of k at each face.
Evaluating k_f*snGrad(T) in my modified solver as you describe:
surfaceScalarField k_f=fvc::interpolate(k);
surfaceScalarField snGradT_f=fvc::snGrad(T);
F=fvc::interpolate(k)*fvc::snGrad(T);
returns the expected fluxes without wiggles. Thank you for this solution.

Can I use OpenFOAMs flux() functions to get the flux used in the laplacian(k,T) scheme immediately, like fvc::flux(phi,T) gets the convective flux?
This is obviously nonsense:
fvc::flux(fvc::laplacian(k, T));
but does a similar function exist, which always returns the used flux independent of the chosen discretisation scheme?

Thanks
Daniel
danielr is offline   Reply With Quote

Old   October 5, 2009, 16:05
Default
  #4
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
Hi Daniel,

I don't think there's a flux function for the diffusive fluxes, at least I didn't find one in the Doxygen djungle (I looked through the fvc namespace).

I note that the flux F you compute above (F=fvc::interpolate(k)*fvc::snGrad(T)) is a surfaceScalarField. Are you ok with that? If you really want something in cell centres to look at in the post-pro, then you can "reconstruct" a volVectorField (!) from your face field F like this:
Fc = fvc::reconstruct(F*mesh.magSf()),
where the surfaceScalarField mesh.magSf() is the face areas.

If you make sure that you use exactly the same fvSchemes for interpolation and snGrad in both the laplacian and the explicit calculation of F, then this Fc will be a very nice vectorfield, even on unstructured grids! Probably the best you can ever hope for. It's in some way "consistent" with your conservative face fluxes, although the notion of "conservation" only makes sense when you look at fluxes in and out of a control volume (like a finite volume cell!), and not vector values in points...

Happy hacking!

/Ola
olwi 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
Wind turbine simulation Saturn CFX 34 October 16, 2014 05:27
Heterogeneous heat flux Kev STAR-CD 4 July 7, 2009 04:48
UDS flux and boundaries Kasper Skriver FLUENT 0 April 18, 2006 08:36
Definition of Mass flux F1 Gavin CD-adapco 1 April 17, 2006 17:51
Replace periodic by inlet-outlet pair lego CFX 3 November 5, 2002 21:09


All times are GMT -4. The time now is 23:43.