the difference between div(U) and div(phi)

kagen816 · July 29, 2022, 05:44

Hi Former,

I am very confused about the difference between div(U) and div(phi). I checked them in a flow field. The difference is about 10^(4to5 to) times. I don't understand why there is a difference. I read some posts and I understand div(phi) is calculating the surface integration of flux, which should be zero ideally. But then, what is the physical meaning of div(U)? shouldn't they be equal to each other? div(U) = div(phi)?

dlahaye · July 29, 2022, 06:08

phi = rho*U

kagen816 · July 29, 2022, 08:37

Quote:

Originally Posted by dlahaye

phi = rho*U

Thanks.Domenico,

But it's not the reason why the difference between div(U) and div(phi) is so big. For the incompressible flow, OpenFOAM rho= 1. I'm still confused.

dlahaye · July 29, 2022, 12:52

U is evaluated on cell centers (cell centered discretization)

phi is evaluated on face center (being a flux).

please check dimensions of U and phi.

kagen816 · July 29, 2022, 13:46

Quote:

Originally Posted by dlahaye

U is evaluated on cell centers (cell centered discretization)

phi is evaluated on face center (being a flux).

please check dimensions of U and phi.

Hi Domenico,

phi's dimension is [0 3 -1 0 0 0 0] (m^3/s)
U's is [0 1 -1 0 0 0 0] (m/s)
But what strange is the dimensions of div(U) and div(phi) is the same. [0 0 -1 0 0 0 0] (/s).

Bloerb · July 30, 2022, 04:21

compressible solvers: phi=rho*U*A
incompressible solvers (since p is defined as p/rho for better performance): phi=U*A
phi is the flux on the surfaces between cells. Hence the unit is either kg/s or m³/s

a div operation basically changes the unit by 1/m. Hence the unit of div(U) is 1/s. The divergence is however defined as follows:
$\nabla\cdot\mathbf{U}=\lim_{V\to 0} \frac{\text{phi}}{|V|}$

It's the sum of the flux over all faces of a cell, as that cell shrinks to zero volume. That should explain the rest

phi=m³/s divided by the volume m³ gives you 1/s as well...

kagen816 · August 2, 2022, 04:32

Quote:

Originally Posted by Bloerb

compressible solvers: phi=rho*U*A
incompressible solvers (since p is defined as p/rho for better performance): phi=U*A
phi is the flux on the surfaces between cells. Hence the unit is either kg/s or m³/s

a div operation basically changes the unit by 1/m. Hence the unit of div(U) is 1/s. The divergence is however defined as follows:
$\nabla\cdot\mathbf{U}=\lim_{V\to 0} \frac{\text{phi}}{|V|}$

It's the sum of the flux over all faces of a cell, as that cell shrinks to zero volume. That should explain the rest

phi=m³/s divided by the volume m³ gives you 1/s as well...

Hi Bloerb,

It's clear for the explanation of the unit. Thank you.

Then, it seems the div operator in OF has more than one definition.

div(U)=∂u/∂x+∂v/∂y+∂w/∂z (U=(u,v,w)) (1)
div(phi)=surface_sum(phi)/cell volume (2)

Then it's hard to understand why there is a difference between the surface_sum(phi) and div(U) because according to Gauss's theorem, they should be the same, shouldn't they?

for (1), I have verified it in Openfoam by calculating the tensor grad(U). div(U) equals grad(U)11+grad(U)22+grad(U)33. it has a big difference with div(phi).

So these two variables have the same unit but different values. is it because the U in the cell centres is not fine enough to get an accurate enough gradient vector (∂u/∂x,∂v/∂y,∂w/∂z)?

Bloerb · August 3, 2022, 01:37

No. These two definitions are identical, check a basic math text book, or wikipedia. Just read up on the definition of these operators. You need to consider, that phi is not a vector field. The divergence is defined for a vector field. It's defintion with the derivatives is based on a vector field. so div(vector) = scalar...which is what i wrote above. div(u) = phi/V. You can express div(U) as your derivatives statement and the right hand side does not contain a div(phi).

Phi is defined as a flux field. A surface scalar field. It contains an area and is a scalar. And you would first need to interpolate a vector field from that and calculate the divergence afterwards. OR you use the definition of the divergence to see directly what that contains. In other words div(phi) itself is not really defined. It takes extra steps to let it make sense.

Hence your difference. You are missing the area somewhere.

kagen816 · August 4, 2022, 03:53

Quote:

Originally Posted by Bloerb

No. These two definitions are identical, check a basic math text book, or wikipedia. Just read up on the definition of these operators. You need to consider, that phi is not a vector field. The divergence is defined for a vector field. It's defintion with the derivatives is based on a vector field. so div(vector) = scalar...which is what i wrote above. div(u) = phi/V. You can express div(U) as your derivatives statement and the right hand side does not contain a div(phi).

Phi is defined as a flux field. A surface scalar field. It contains an area and is a scalar. And you would first need to interpolate a vector field from that and calculate the divergence afterwards. OR you use the definition of the divergence to see directly what that contains. In other words div(phi) itself is not really defined. It takes extra steps to let it make sense.

Hence your difference. You are missing the area somewhere.

Hi Bloerb,

Thank you.

I have checked the div(U) by a very simple model. i put some monitors at the centres of the cells and outputed the velocity. I interpolated the cell centre data to the cell surface centres. and then, I picked one cell to calculate the flux sum of cell surfaces and divided by cell volume: sum(Uf&Sf)/cell_volume.

The results give a number like 0.125 and it's exactly the same as the direct output fvc::div(U) from OF. So now I know how OF calculate the fvc::div(U).

However, for the fvc::div(phi), the output value from OF is like -4.7e-8, which is much smaller than fvc::div(U). and I totally have no idea how this number is calculated.

In your answer above you mentioned '' you would first need to interpolate a vector field from that and calculate the divergence afterwards'', do you mean I should create a vector field from the scalar field?

Bloerb · August 4, 2022, 12:14

OK, a bit more basics because i was a bit sloppy in my explanations....As you stated this is how the divergence is calculated:

$\nabla\cdot \mathbf{u} = \frac{\partial u_x}{\partial x}+ \frac{\partial u_y}{\partial y}+ \frac{\partial u_z}{\partial z}$

There is however a BIG BIG problem with that...First of all..this is for continous fields...we are however using discretized cells...hence none of these derivatives are calculated without additional work. You are only saving the velocity at the center of each cell, and the flux field phi on the surface of each cell.....now as i stated before, the formal definition of the divergence stems from vector calculus. The divergence states, that the integral of the divergence of a vector field over a certain volume can simply be stated as the integral of the same thing over the surface. The divergence does hence reduce a volume integral to a surface integral...this is called the divergence theorem or Gauss theorem and is taught in every fluid dynamics course. It is basically this:

$\nabla\cdot \mathbf{u} = \frac{1}{V}\int_V\left(\nabla\cdot \mathbf{u} \right)dV = \frac{1}{V}\oint_S\left(\mathbf{n}\cdot\mathbf{u}\right)dS = \frac{1}{V}\sum_i^{nFace}u_{f,i}\cdot s_{f,i}$

Now the beauty of this is that instead of calculating the divergence by calculating the derivatives you simply sum up all the fluxes phi=u_face*A_face*n=u*s....and divide by the cells volume...think about that for a second...on an arbitrary tet or polyeder mesh for example...how do you calculate du/dx? Not that simple....On a normal hex mesh you'd need to take the velocity at two cells centers and calculate the difference in x direction between them and divde by the distance...but which cell the one to the left or the one on the right?...and which cell do you use on a tet or poly mesh? none of the cells are oriented in x direction....which cell do you use...you potentially need more than one.... Not that simple....and leads to lots of numerical error potential...hence OpenFOAM uses the Gauss theorem for calculating the divergence. Just sum up phi and divide by volume. Which reduces the complexity enormously. This is also in line with the use case for divergence in fluid dynamics...the mass balance equation is div(U)=0 for incompressible media...and this is as the equation shows...nothing but, the sum of every flow entering/leaving the cell...and this needs to be zero.

As is stated before....divergence of a scalar is mathematically not defined...Hence div(phi) and div(U) are not identical, nor should they...div(phi) itself does not exist!

Now if you want to get into the details...why div(phi) exists in OpenFOAM check fvcDiv.C and fvcSurfaceIntegrate divScheme and convectionScheme.C/H. Please keep in mind that openFOAMs div(phi,U) etc are methods for calculating the volume integral of an equation... div(u) is a method for calculating 1/V int div(u) dv ... not for div(u)...and
div(phi) uses surface integrate...hence the difference you are experiencing should be the area.

fly_light · August 5, 2023, 15:42

Hello Kagen,
I have the same problem with the different results of div(U) and div(phi)

Quote:

However, for the fvc::div(phi), the output value from OF is like -4.7e-8, which is much smaller than fvc::div(U)

Did you find the answer?
I read the source code of fvc::div class, and the value is calculated by fvc::surfaceIntegrate, but what I understand from the code is that the two values are the same. I must make some mistakes.
Please tell me if you find the answer.

Best regards.

July 29, 2022, 05:44	the difference between div(U) and div(phi)	#1
kagen816 New Member Not applicable Join Date: Sep 2018 Posts: 15 Rep Power: 8	Hi Former, I am very confused about the difference between div(U) and div(phi). I checked them in a flow field. The difference is about 10^(4to5 to) times. I don't understand why there is a difference. I read some posts and I understand div(phi) is calculating the surface integration of flux, which should be zero ideally. But then, what is the physical meaning of div(U)? shouldn't they be equal to each other? div(U) = div(phi)?

July 30, 2022, 04:21		#6
Bloerb Senior Member Join Date: Sep 2013 Posts: 353 Rep Power: 21	compressible solvers: phi=rhoUA incompressible solvers (since p is defined as p/rho for better performance): phi=U*A phi is the flux on the surfaces between cells. Hence the unit is either kg/s or m³/s a div operation basically changes the unit by 1/m. Hence the unit of div(U) is 1/s. The divergence is however defined as follows: $\nabla\cdot\mathbf{U}=\lim_{V\to 0} \frac{\text{phi}}{\|V\|}$ It's the sum of the flux over all faces of a cell, as that cell shrinks to zero volume. That should explain the rest phi=m³/s divided by the volume m³ gives you 1/s as well... jherb, hogsonik, kagen816 and 1 others like this.

August 3, 2022, 01:37		#8
Bloerb Senior Member Join Date: Sep 2013 Posts: 353 Rep Power: 21	No. These two definitions are identical, check a basic math text book, or wikipedia. Just read up on the definition of these operators. You need to consider, that phi is not a vector field. The divergence is defined for a vector field. It's defintion with the derivatives is based on a vector field. so div(vector) = scalar...which is what i wrote above. div(u) = phi/V. You can express div(U) as your derivatives statement and the right hand side does not contain a div(phi). Phi is defined as a flux field. A surface scalar field. It contains an area and is a scalar. And you would first need to interpolate a vector field from that and calculate the divergence afterwards. OR you use the definition of the divergence to see directly what that contains. In other words div(phi) itself is not really defined. It takes extra steps to let it make sense. Hence your difference. You are missing the area somewhere. Last edited by Bloerb; August 3, 2022 at 16:30.

August 4, 2022, 12:14		#10
Bloerb Senior Member Join Date: Sep 2013 Posts: 353 Rep Power: 21	OK, a bit more basics because i was a bit sloppy in my explanations....As you stated this is how the divergence is calculated: $\nabla\cdot \mathbf{u} = \frac{\partial u_x}{\partial x}+ \frac{\partial u_y}{\partial y}+ \frac{\partial u_z}{\partial z}$ There is however a BIG BIG problem with that...First of all..this is for continous fields...we are however using discretized cells...hence none of these derivatives are calculated without additional work. You are only saving the velocity at the center of each cell, and the flux field phi on the surface of each cell.....now as i stated before, the formal definition of the divergence stems from vector calculus. The divergence states, that the integral of the divergence of a vector field over a certain volume can simply be stated as the integral of the same thing over the surface. The divergence does hence reduce a volume integral to a surface integral...this is called the divergence theorem or Gauss theorem and is taught in every fluid dynamics course. It is basically this: $\nabla\cdot \mathbf{u} = \frac{1}{V}\int_V\left(\nabla\cdot \mathbf{u} \right)dV = \frac{1}{V}\oint_S\left(\mathbf{n}\cdot\mathbf{u}\right)dS = \frac{1}{V}\sum_i^{nFace}u_{f,i}\cdot s_{f,i}$ Now the beauty of this is that instead of calculating the divergence by calculating the derivatives you simply sum up all the fluxes phi=u_faceA_facen=u*s....and divide by the cells volume...think about that for a second...on an arbitrary tet or polyeder mesh for example...how do you calculate du/dx? Not that simple....On a normal hex mesh you'd need to take the velocity at two cells centers and calculate the difference in x direction between them and divde by the distance...but which cell the one to the left or the one on the right?...and which cell do you use on a tet or poly mesh? none of the cells are oriented in x direction....which cell do you use...you potentially need more than one.... Not that simple....and leads to lots of numerical error potential...hence OpenFOAM uses the Gauss theorem for calculating the divergence. Just sum up phi and divide by volume. Which reduces the complexity enormously. This is also in line with the use case for divergence in fluid dynamics...the mass balance equation is div(U)=0 for incompressible media...and this is as the equation shows...nothing but, the sum of every flow entering/leaving the cell...and this needs to be zero. As is stated before....divergence of a scalar is mathematically not defined...Hence div(phi) and div(U) are not identical, nor should they...div(phi) itself does not exist! Now if you want to get into the details...why div(phi) exists in OpenFOAM check fvcDiv.C and fvcSurfaceIntegrate divScheme and convectionScheme.C/H. Please keep in mind that openFOAMs div(phi,U) etc are methods for calculating the volume integral of an equation... div(u) is a method for calculating 1/V int div(u) dv ... not for div(u)...and div(phi) uses surface integrate...hence the difference you are experiencing should be the area. jherb, ashvinc9, acgnipper and 2 others like this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
mass imbalance / div(U) vs div(phi)	geth03	OpenFOAM Post-Processing	3	August 5, 2023 15:14
divergence evaluation: direct (div(U)) or through face fluxes (div(phi))	t.teschner	OpenFOAM Programming & Development	0	March 16, 2018 13:26
difference between div(U) and div(phi)	fshak92	OpenFOAM Programming & Development	2	August 15, 2014 08:53

July 29, 2022, 06:08		#2
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 772 Blog Entries: 1 Rep Power: 17	phi = rho*U

July 29, 2022, 12:52		#4
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 772 Blog Entries: 1 Rep Power: 17	U is evaluated on cell centers (cell centered discretization) phi is evaluated on face center (being a flux). please check dimensions of U and phi.