CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   discrepancy in fvc::grad(U) (https://www.cfd-online.com/Forums/openfoam-programming-development/219652-discrepancy-fvc-grad-u.html)

backscatter August 3, 2019 19:26

discrepancy in fvc::grad(U)
 
HI all,

Recently I noticed that in openfoam, gradient of a vector is computed as the outer product of the gradient operator and the vector. Page 25 (bottom) of the programmer's guide defines the gradient of velocity field as:
Code:

\nabla_{i}a_{j}
which is transpose of what we exactly mean by grad(U) as:
Code:

\nabla_{j}a_{i}
I checked the computed values of gradient of velocity field near the wall (where wall normal component is expected to be dominant) and I can confirm that it outputs the transpose of what we actually want.
To my surprise, I have seen people use the grad(U) without taking its transpose in openfoam that would actually give the correct tensor.
Has anybody encountered this issue before?
OpenFOAM Version: 2.4.0

anon_q August 3, 2019 21:11

Hello,
The result is correct however the notation of the manual is confusing (wrong).
The problem is due to the nabla operator (\nabla).
because by convention [Ref. 1] (look carefuly to the position of the indices):

\dfrac{\partial a_i}{\partial x_j} \equiv \partial_ja_i
Why? because in this context, a_i\partial_j will not make sense (remember that the operator nabla must act on the vector \textbf{a} so it must be at the left).

Conclusion:
The gradient is calculated correctly (matrix form), but the notation (\partial_i a_j) given in the manual for the gradient is wrong (it is for the transpose of the gradient).

[Ref. 1]: See page 7 of this document: http://web.iitd.ac.in/~pmvs/courses/mcl702/notation.pdf

backscatter August 3, 2019 21:31

Quote:

Originally Posted by Evren Linda (Post 740999)
Hello,
The result is correct however the notation is confusing.
The problem is due to the nabla operator (\nabla).
because by convention (look carefuly to the position of the indices):

\dfrac{\partial a_i}{\partial x_j} \equiv \partial_ja_i
Why? because in this context, a_i\partial_j will not make sense (remember that the operator nabla must act on the vector \textbf{a} so it must be at the left).

Conclusion:
The gradient is calculated correctly, but the notation given in the manual is wrong (it is for the transpose of the gradient).

I checked the resulting tensor in a simple channel flow where dU/dy is the dominant term. So we would expect the second term (XY) to be large compared to others. This is what I found where it is very clear that the grad(U) computed by the fvc:grad() function is actually the transpose of the actual gradient.
The computed velocity gradient tensor computed by grad() operator for channel flow (near wall) is below:
(0.017149094 -0.00010358519 -0.016898076 29.161076 0.0014569992 1.2603341 0.040551971 -0.00023234273 0.013089031)
As you can see, the fourth component (YX) in red is very clearly the wall normal component as it is the dominant term which should actually be the second (XY) component.
I guess whatever is presented in the guide is correct.
So if the gradient computed is correct, shouldn't the dominant (dU/dy) term be second in the above tensor?

anon_q August 3, 2019 21:41

The result is correct because the gradient will be calculated as follows:

\begin{bmatrix}
\partial_1 U & \partial_1 V & \partial_1 W \\
{\color{green}\partial_2 U} & \partial_2 V & \partial_2 W \\
\partial_3 U & \partial_3 V & \partial_3 W 
\end{bmatrix}

As you can see: {\color{green}\partial_2 U} \equiv \dfrac{\partial U}{\partial y} is the 4th element in the tensor (correctly calculated by OpenFOAM).
Remember that the tensor class will represent it like this:
gradU =Tensor(\partial_1 U, \partial_1 V, \partial_1 W, {\color{green}\partial_2 U}, \partial_2 V, \partial_2 W, \partial_3 U, \partial_3 V, \partial_3 W) but to access {\color{green}\partial_2 U \equiv \dfrac{\partial U}{\partial y}} we use gradU.yx()

backscatter August 3, 2019 21:55

Quote:

Originally Posted by Evren Linda (Post 741001)
The result is correct because the gradient will be calculated as follows:

\begin{bmatrix}
\partial_1 U & \partial_1 V & \partial_1 W \\
{\color{green}\partial_2 U} & \partial_2 V & \partial_2 W \\
\partial_3 U & \partial_3 V & \partial_3 W 
\end{bmatrix}

As you can see: {\color{green}\partial_2 U} \equiv \dfrac{\partial U}{\partial y} is the 4th element in the tensor (correctly given by OpenFOAM)

Yes, absolutely. That's what my point is. The programmers guide is correct and unlike common definition of gradU in textbooks, the one that openfoam computes is its transpose.
This is what we actually need (the common agreed upon definition of velocity gradient in all textbooks on fluid mechanics):
Code:

U1,1 U1,2 U1,3
U2,1 U2,2 U2,3
U3,1 U3,2 U3,3

And this is what openfoam computes:
Code:

U1,1 U2,1 U3,1
U1,2 U2,2 U3,2
U1,3 U2,3 U3,3

which is clearly the transpose of what we usually mean by gradient of any vector field in textbooks.

anon_q August 3, 2019 22:21

This is a matter of layout convention. You can take a look to this detailed section on wikipedia.

Your notation:
Quote:

U1,1 U1,2 U1,3
U2,1 U2,2 U2,3
U3,1 U3,2 U3,3
is commonly used for the jacobian which is the transpose of the gradient.

backscatter August 3, 2019 22:31

Quote:

Originally Posted by Evren Linda (Post 741004)
This is a matter of layout convention. You can take a look to this detailed section on wikipedia.

Your notation:


is commonly used for the jacobian which is the transpose of the gradient.

You're right.
The issue was just to make it very known that with using fvc::grad(U) in foam, as I said, we are not computing the gradient as the textbook definition used for gradient of velocity is written.
I have seen people use it right away in calculations when they actually were meant to use its transpose. So care must be taken!
Thanks for your time!! I am pretty sure that this discussion will help a lot of people!

anon_q August 3, 2019 23:07

In my opinion, I find the matrix form of the gradient to be correct. however the notation \partial_ia_j for the gradient is wrong.

backscatter August 3, 2019 23:25

Quote:

Originally Posted by Evren Linda (Post 741007)
In my opinion, I find the matrix form of the gradient to be correct. however the notation \partial_ia_j for the gradient is wrong.

No it is not wrong. It expands to a tensor that you mentioned in your second comment to this thread which is the correct form of gradient.

fanny June 8, 2021 04:35

Hello,

I also encountered this issue. In my opinion, fvc::grad computes the transpose of the gradient... It is well documented in the programmers' guide - p25 eq2.3.

Tobermory June 10, 2021 13:05

No, I think that the OpenFOAM notation does hold up, and makes more sense to me.

For example, the outer product of two vectors, \mathbf{P\mathbf} = \mathbf{a} \otimes \mathbf{b}, can be written in terms of matrix multiplication as \mathbf{P} = \mathbf{a} \mathbf{b}^T.
Now apply that to the grad operation, which is the outer product of nabla and the vector:

\mathbf{G} = grad \:\mathbf{b} = \nabla \otimes \mathbf{b}
=
 \left( \begin{array}{c} \partial/\partial x \\ \partial/\partial y \\ \partial/\partial z \end{array} \right) \otimes
 \left( \begin{array}{c} b_x \\b_y \\ b_z \end{array} \right)
=
 \left( \begin{array}{c} \partial/\partial x \\ \partial/\partial y \\ \partial/\partial z \end{array} \right) 
 \left( \begin{array}{ccc} b_x & b_y & b_z \end{array} \right)

which on multiplying out gives:

\mathbf{G} =  \left( \begin{array}{ccc}
 \partial b_x/\partial x & \partial b_y/\partial x &  \partial b_z/\partial x \\ 
 \partial b_x/\partial y & \partial b_y/\partial y &  \partial b_z/\partial y \\ 
 \partial b_x/\partial z & \partial b_y/\partial z &  \partial b_z/\partial z \\ 
 \end{array} \right)

and hey presto you get the convention that OpenFOAM uses. Just as a reminder, the middle term on the top row is component \mathbf{G}_{12}= \partial /\partial x (b_y) = \partial b_y/\partial x. You just need to remember that the nabla operator comes before the thing it operates on, and so its index comes first.

backscatter November 22, 2021 23:38

Also, we are not much bothered about the convention. We are more concerned about the correct implementation of the terms that arise in FD equations. I think the best practice to address this is to be careful what the term that you are trying to implement involves -- whether you want [U1,1 U1,2 U1,3; U2,1 U2,2 U2,3; U3,1 U3,2 U3,3 ] to be implemented or its transpose. The implementation of gradient in FOAM is technically its transpose but in most of the fluid dynamic equations the tensor I just wrote above is more common representation of (del/del_xj U_i).

backscatter November 22, 2021 23:47

As an example, Tobermory and Fanny, how would you implement the following term:

\overline{u_i u_k} \frac{\partial U_j}{\partial x_k} + \overline{u_j u_k} \frac{\partial U_i}{\partial x_k}

Please submit your responses below in terms of expanded tensors. Thanks

Tobermory November 24, 2021 18:06

1 Attachment(s)
Quote:

Originally Posted by backscatter (Post 817115)
Also, we are not much bothered about the convention. We are more concerned about the correct implementation of the terms that arise in FD equations.

In my mind, there is no "correct" way to write it ... just convention. And opinions seem to be split on this. Some sources write it the way you do, others write it the way that OpenFOAM writes it, i.e. the transpose (see eg page 7 of https://web.iitd.ac.in/~pmvs/courses...2/notation.pdf). I have explained in my earlier post why I think the OpenFOAM way makes sense, to me. But so long as you use the tensor in a consistent fashion, then it does not matter, I think.

As for your Reynolds stress production terms, see the attached image for the OpenFOAM notation way of writing this (I ran out of patience with trying to do the notation in this editor!).

Hope that helps ...

Note: some of the terms in the last matrix are wrong - see the corrected version below

backscatter November 24, 2021 18:25

Quote:

Originally Posted by Tobermory (Post 817313)
In my mind, there is no "correct" way to write it ... just convention. And opinions seem to be split on this. Some sources write it the way you do, others write it the way that OpenFOAM writes it, i.e. the transpose (see eg page 7 of https://web.iitd.ac.in/~pmvs/courses...2/notation.pdf). I have explained in my earlier post why I think the OpenFOAM way makes sense, to me. But so long as you use the tensor in a consistent fashion, then it does not matter, I think.

As for your Reynolds stress production terms, see the attached image for the OpenFOAM notation way of writing this (I ran out of patience with trying to do the notation in this editor!).

Hope that helps ...

Yes, you are right about there being no "correct" way to write it. It is just about getting the implementation right.
Correct, that is how I would implement the first term. Would the second term implementation be any different?
I would write the second term as: (with the velcoity gradient getting pre-multiplied)
\frac{\partial U_i}{\partial x_k} \overline{u_j u_k}
with the stress tensor and your velocity gradient tensor transposed as:

\begin{pmatrix}U1,1&U1,2&U1,3\\ U2,1&U2,2&U2,3\\ U3,1&U3,2&U3,3\end{pmatrix}\begin{pmatrix}uu&vu&wu\\ uv&vv&wv\\ uw&vw&ww\end{pmatrix}

Do you agree?

Tobermory November 25, 2021 09:37

1 Attachment(s)
Yes - exactly right. The second Re stress production term \overline{u_j u_k} \frac{\partial U_i}{\partial x_k} is just the transpose of the first - check out the attached figure (I also corrected the first term from my earlier post).

Voulet November 27, 2021 08:01

According to my own test. I confirm that OpenFoam uses the transpose of the jacobian matrix for the gradient of a vector.


I think that it should be written somewhere because when you derive your own equation using the classical jacobian matrix it can be a source of error very hard to catch if you code it without knowing it.

HPE March 23, 2022 04:17

The treatment on the grad is confusing in my opinion.

Just imagine you compute "skew(gradU)" by assuming the traditional notation.


All times are GMT -4. The time now is 20:30.