CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Harmonic average of Diffusion Tensors in Finite Volume Method

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 5, 2018, 04:26
Default Harmonic average of Diffusion Tensors in Finite Volume Method
  #1
New Member
 
Join Date: Sep 2018
Posts: 3
Rep Power: 4
Yurie_Breschnef is on a distinguished road
Hi Everyone, this is my first post on this site so please be gentle.

I want to implement a Bilinear Finite Volume discretization
of the anisotropic diffusion problem:


\frac{du}{dt} = \nabla \cdot (\textbf{D} \nabla u)

Both my degrees of freedom, as well as the Diffusion Tensors, are given
on the vertices of my quadrilateral grid. I can, therefore, use bilinear
test/trialfunctions on each cell to calculate the full gradient.

I am wondering and am deeply puzzled as to how I am to approximate the
Diffusion tensor on an internal sub-face of my control volume.

In the 1D case (TPFA), the diffusion coefficient is modeled to be constant
over one cell and the degree of freedom is thought to be the same
everywhere within the cell. The flux over an interface is made consistent
by using the harmonic average of the two given diffusivities:


f_{A->B}= \frac{2}{\frac{1}{D_A}+ \frac{1}{D_B}} ~\nabla_x u

Usually, the gradient is then approximated by the simple finite difference
of the values modeled to be in the cell centers.

sketch to clarify:



So far so good. Now I am not in a setting where the Grid is aligned with
the principal directions of the diffusion tensors in involved (not "K-orthogonal").
The Two-point flux approximation would, therefore, fail, because it can not
capture the off-diagonal entries of the diffusion tensors involved, because
only the x-derivative can be calculated.

In my 2D setting, **I do reconstruct the full gradient**. I want to calculate
four sub-fluxes within my cell as indicated by the following sketch:




Here the black dots indicate where my degrees of freedom lie. The dashed line shows the area where I model the diffusion tensor to be constant (dual grid), the friendly green dot marks the location where I want to evaluate the bottom-most of the four sub-fluxes to properly assemble. Following the procedure from the 1D case, I could use the harmonic average of the diffusion tensors of the involved dual cells A,B. Let's call it

**option 1:**

D_{eff} = \big(\frac{\textbf{A}^{-1 } + \textbf{B}^{-1} }{2}\big)^{-1}

Here the "-1" in the exponents indicates the matrix inverse. What puzzles me so much is that by calculating the matrix inverse, the diagonal entries become coupled with the off-diagonal ones. To my understanding, these are different physical processes and I feel awkward to have, say, The component D_{xx} depend on D_{xy} etc.

The second option I see is to do a harmonic average component-wise, to have:

**option 2**

D_{eff} = \begin{bmatrix}
 \frac{2}{\frac{1}{D^{A}_{xx}} + \frac{1}{D^{B}_{xx}}} &
\frac{2}{\frac{1}{D^{A}_{xy}} + \frac{1}{D^{B}_{xy}}} \\
\frac{2}{\frac{1}{D^{A}_{yx}} + \frac{1}{D^{B}_{yx}}} &
\frac{2}{\frac{1}{D^{A}_{yy}} + \frac{1}{D^{B}_{yy}}} 
 \end{bmatrix}

If i do it this way, the diagonal component D_{xx} is only dependend on the corresponding entries in D_A and D_B.


Investigating option two I found that as the offdiagonal entries of the Diffusion Matrices may be negative. So I face the situation where I take the harmonic average of a negative and a positive number which, according to wikipdeia, is not defined properly. So now both methods do not seem to make sense.


What is the correct way to do average the Diffusion Tensors at the sub-face marked with the green dot?

[1]: https://i.stack.imgur.com/LGl4N.png
[2]: https://i.stack.imgur.com/9BcZQ.png
Yurie_Breschnef is offline   Reply With Quote

Old   September 5, 2018, 11:29
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,687
Rep Power: 60
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
I suppose you have D.Grad(u), therefore a vector. Then you apply the divergence and get a scalar function.

Are all the variable co-located? I don't see specific problems due to the tensor.
You can also think to work with a FV approach starting from the integral form.
FMDenaro is offline   Reply With Quote

Old   September 6, 2018, 04:06
Default
  #3
New Member
 
Join Date: Sep 2018
Posts: 3
Rep Power: 4
Yurie_Breschnef is on a distinguished road
Thank you for answering!
I am in a finite Volume Setting, where I apply the gauss-theorem to write the volume integral over: nabla (D.grad(u)) as a integral over the faces. There ,to evaluate the fluxes over the boundaries I would do the scalar product with the normal of the faces (the indicated subface in the above post). The flux calculated can of course be negative, so again I face the situation of doing the harmonic average over negative values. It seems that I have to find an equivalent of to the harmonic average of the diffusion coefficient in the 1D case, but for the whole matrix (2D).
Yurie_Breschnef is offline   Reply With Quote

Old   September 6, 2018, 04:26
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,687
Rep Power: 60
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Yurie_Breschnef View Post
Thank you for answering!
I am in a finite Volume Setting, where I apply the gauss-theorem to write the volume integral over: nabla (D.grad(u)) as a integral over the faces. There ,to evaluate the fluxes over the boundaries I would do the scalar product with the normal of the faces (the indicated subface in the above post). The flux calculated can of course be negative, so again I face the situation of doing the harmonic average over negative values. It seems that I have to find an equivalent of to the harmonic average of the diffusion coefficient in the 1D case, but for the whole matrix (2D).

Using the integral form you have


d/dt Int[V] u dV = Int[S] n.(D.grad u) dS


therefore when you sum all the fluxes, the RHS must be negative for a diffusive problem. That does not say the sign of the flux on each face. And what you compute in the RHS is always an average (just divide by the volume measure).

Does the matrix entries depend on u?
FMDenaro is offline   Reply With Quote

Old   September 6, 2018, 04:51
Default
  #5
New Member
 
Join Date: Sep 2018
Posts: 3
Rep Power: 4
Yurie_Breschnef is on a distinguished road
Hi!
The Diffusion matrices are given by a dataset (DTI) and are symmetric-positive definite and constant over the simulation. They may have strong inhomogeneities so there might be orders of magnitude difference between A and B. I want the subflux (indicated by the green dot) which I calculate seen from cell A:

F_a = \mathbf{n} \cdot(\mathbf{D_a}  \nabla u)*h

to be consistent with the Flux seen from cell B (bottom right):

F_b = \mathbf{n}\cdot(\mathbf{D_b}  \nabla u)*h
(here h is the length(2D)/area(3D) of the subface)

But D_a might be very different to D_b.
Since I evaluate the gradient by bilinear approximation, it is the same in both calculations of the fluxex F_a and F_b, the normal is also identical. I need a way to 'average' the diffusion tensors before I apply them to the gradient. If I calculate the full fluxes F_a and F_b and then do a harmonic average, the result does not make sense for negative fluxes.
I need something like:

\mathbf{D_{eff}} = AVG(\mathbf{D_a},\mathbf{D_b})

So I can consistently discretise my flux over the subface (for both cells):

F_{consistent} =\mathbf{n}\cdot(\mathbf{D_{eff}}  \nabla u)*h


In the 1D case this average is the mentioned harmonic mean.

the two options I showed do not make sense to me. Thanks for taking the time!
Yurie_Breschnef is offline   Reply With Quote

Reply

Tags
diffusion coefficient, finite-volume, harmonic mean

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lattice Boltzmann method vs Finite Element Method and Finite Volume Method solemnpriest Main CFD Forum 3 August 12, 2013 11:00
Control volume based finite difference method? mukut.medhi Main CFD Forum 3 August 24, 2012 10:01
Finite Volume Method cfd seeker Main CFD Forum 3 September 8, 2011 04:36
Time unit in finite volume solution for diffusion problem. alpha2beta Main CFD Forum 0 June 30, 2010 02:05
fluent add additional zones for the mesh file SSL FLUENT 2 January 26, 2008 11:55


All times are GMT -4. The time now is 15:49.