discrepancy in fvc::grad(U)

backscatter · August 3, 2019, 20:26

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:

which is transpose of what we exactly mean by grad(U) as:

Code:

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, 22: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, 22:31

Quote:

Originally Posted by Evren Linda

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, 22: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, 22:55

Quote:

Originally Posted by Evren Linda

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, 23: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, 23:31

Quote:

Originally Posted by Evren Linda

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 4, 2019, 00: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 4, 2019, 00:25

Quote:

Originally Posted by Evren Linda

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, 05: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, 14: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 23, 2021, 00: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 23, 2021, 00: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, 19:06

Quote:

Originally Posted by backscatter

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, 19:25

Quote:

Originally Posted by Tobermory

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, 10:37

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, 09: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, 05:17

The treatment on the grad is confusing in my opinion.

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

August 3, 2019, 20:26	discrepancy in fvc::grad(U)	#1
backscatter Member Anonymous Join Date: Aug 2016 Posts: 75 Rep Power: 10	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 C. Okubo, Voulet, HPE and 2 others like this.

August 3, 2019, 22:11		#2
anon_q Senior Member Join Date: Mar 2018 Posts: 115 Rep Power: 8	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 C. Okubo likes this.

August 3, 2019, 22:41		#4
anon_q Senior Member Join Date: Mar 2018 Posts: 115 Rep Power: 8	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() watermelon likes this.

August 4, 2019, 00:07		#8
anon_q Senior Member Join Date: Mar 2018 Posts: 115 Rep Power: 8	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. lpz456 likes this.

June 8, 2021, 05:35		#10
fanny New Member Join Date: Feb 2016 Posts: 13 Rep Power: 10	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. lpz456 likes this.

June 10, 2021, 14:05		#11
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 736 Rep Power: 14	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. lpz456 likes this. Last edited by Tobermory; June 14, 2021 at 07:50. Reason: Tidying up the wording.

November 23, 2021, 00:38		#12
backscatter Member Anonymous Join Date: Aug 2016 Posts: 75 Rep Power: 10	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). lpz456 likes this.

November 23, 2021, 00:47		#13
backscatter Member Anonymous Join Date: Aug 2016 Posts: 75 Rep Power: 10	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

November 27, 2021, 09:01		#17
Voulet Member Join Date: Jun 2019 Posts: 41 Rep Power: 7	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. saladbowl, HPE, openfoam_aero and 1 others like this. __________________ « Debugging is what CFD is about. 5 minutes to modify your code, 5 months to find why it does not work anymore. »

March 23, 2022, 05:17		#18
HPE Senior Member Herpes Free Engineer Join Date: Sep 2019 Location: The Home Under The Ground with the Lost Boys Posts: 931 Rep Power: 13	The treatment on the grad is confusing in my opinion. Just imagine you compute "skew(gradU)" by assuming the traditional notation. __________________ The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki-1/Wiki-2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESI-OpenCFD-Trademark) Bugs/Features: FOAM-Extend (Wikki-FSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror Last edited by HPE; March 23, 2022 at 07:37.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
how can I do dot product of U and fvc::grad(U) where U is an volVectorFiled ?	mechy	OpenFOAM	3	February 24, 2023 12:17
Acoustic Power Discrepancy	Time4Tea	FLUENT	0	August 21, 2018 10:46
implicit - (fvc::grad(U) & fvc::grad(muEff))	Kareem Abdelshafy	OpenFOAM Programming & Development	1	August 20, 2016 21:57
strange behaviour of fvc::grad(U) and NonlinearKEShih turbmod - analitical comparison	ancsa	OpenFOAM Verification & Validation	1	January 20, 2015 12:00
FLUENT Version Mis-Match Discrepancy	asd	Main CFD Forum	0	April 19, 2007 05:58