Harmonic average of Diffusion Tensors in Finite Volume Method

Yurie_Breschnef · September 5, 2018, 04:26

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

and

.

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

FMDenaro · September 5, 2018, 11:29

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.

Yurie_Breschnef · September 6, 2018, 04:06

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).

FMDenaro · September 6, 2018, 04:26

Quote:

Originally Posted by Yurie_Breschnef

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?

Yurie_Breschnef · September 6, 2018, 04:51

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!

September 5, 2018, 04:26	Harmonic average of Diffusion Tensors in Finite Volume Method	#1
Yurie_Breschnef New Member Join Date: Sep 2018 Posts: 3 Rep Power: 7	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 and . 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

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

September 5, 2018, 11:29		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	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.

September 6, 2018, 04:06		#3
Yurie_Breschnef New Member Join Date: Sep 2018 Posts: 3 Rep Power: 7	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).

September 6, 2018, 04:51		#5
Yurie_Breschnef New Member Join Date: Sep 2018 Posts: 3 Rep Power: 7	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!