View Issue Details

IDProjectCategoryView StatusLast Update
0003619OpenFOAMBugpublic2021-02-10 12:15
Reportertiziano.maffeiAssigned Tohenry 
Status closedResolutionno change required 
PlatformGNU/LinuxOSUbuntuOS Version18.04
Product Version8 
Fixed in Version 
Summary0003619: Issue using fvOptions with semiImplicitSource option
DescriptionGood morning,

I have carried out a few mixing simulations of two streams with different inlet temperatures and with an energy source placed in a cellZone of the domain.
In order to take into account the energy source, the semiImplicitSource option (based on enthalpy) has been set up in the system/fvOptions file.
The simulations have been performed with several versions of OpenFoam: 5, 6, 7, 8.

At the outlet patch, the temperature weighted on the mass flow rate has been calculated.

The mass flow rate weighted temperature of simulations 5 and 6 is very close to the theoretical temperature. On the other side, even if the energy balance is properly closed, the weighted outlet temperature of simulations 7,8 is different from the theoretical ones. In this case, the difference between 5/6 and 7/8 is around 3°. In more complex geometries, the difference is much higher (30°C).

I haven't found any difference between the versions when the semiImplicitSource option is off.


Steps To ReproduceCase OF6 and OF8 are attached
Each case contains postProcess folder with weighted temperature

To run the case:
- do not delete polyMesh folder (mesh build with Ansys and not attached here due to size limitations)
- solver: rhoSimpleFoam
TagsNo tags attached.



2021-02-03 15:00

updater   ~0011831

@tiziano.maffei: You can attach the packages below after creating the report.
I've closed the other two duplicates for now and will be erased later on.

Notice that the limit per file is 2MB. Please keep in mind to provide a simple-enough case to diagnose, hence the file size limit.


2021-02-03 15:06

reporter   ~0011832

sorry for the inconvenience!

To run the case:
./Allrun ( I have also attached the file *.msh)


TjunctionEnthalpySourceOF8.7z (843,122 bytes)


2021-02-04 14:17

manager   ~0011833

I ran the case in OpenFOAM-dev and it converges and the results look fine. From the information you have given above it is not clear what result you are looking for.

Also the thermo specification is complex, it would make sense to run with much simpler thermo to diagnose the problem, constant Cp, constant transport, constant rho etc.


2021-02-04 14:55

reporter   ~0011834

Dear Henry,
thank you for your answer.

Yes, you right regarding apply simpler thermo properties.
In the beginning, I have performed a few simulations with a simpler thermoDict (with constant Cp, constant transport and perfect gas), but I didn't find any difference between versions 6 and 8. The discrepancy comes up when I select polynomial properties.

My point is why OF6 and OF8 predict a different outlet temperature. In this case, when I run this simulation with OpenFoam 5 or 6, the weighted outlet temperature (based on mass flow rate) is around 614.13 K. If I calculate the outlet temperature by hand, simply making an energy balance, I obtained 614.14 K.
When I run this case with the version OF8, the weighted outlet temperature is around 611.42 K

I hope you can find my point more clear. Sorry for my unclear explanation.



2021-02-04 15:04

manager   ~0011835

If you set Cp as a function of temperature you need to integrate h rather than T to check conservation, the integral of T is not a conserved property.

The difference between version 6 and 8 is in the handling of the departure function. In OpenFOAM-6 the departure function is not implemented for icoPolynomial but it is in OpenFOAM-8 so that now the results from solving for enthalpy and internal energy are consistent. If you are running with an incompressible equation of state it make a lot more sense to solve for internal energy rather than enthalpy.


2021-02-04 15:41

reporter   ~0011836

Dear Henry,
when I have calculated by hand the temperature I have properly calculated the enthalpy, as integral of Cp(T). To find the outlet temperature, I have changed the outlet Temperature automatically in order to guarantee Hin-Hout=0.

I have also performed a simulation with selecting sensible internal energy in thermoDict and internal energy source in system/fvOption with OF8. The result of this simulation is basically the same as the simulation carried out with enthalpy. I have gotten 611.4K but not the "theoretical temperature" of 614.14 K



2021-02-04 15:48

manager   ~0011837

It would make sense to perform integrals on h rather than to check conservation.


2021-02-04 15:54

manager   ~0011838

> "theoretical temperature" of 614.14 K

Where does this "theoretical temperature" value come from? Does the case have an analytical solution?


2021-02-04 16:02

reporter   ~0011839

Dear Henry,

I have improperly called the temperature calculated by hand "theroretical temperature". I have also checked this outlet temperature with a process simulator and I have obtained 614.4K
But is it not enough to check the enthalpy conservation?


2021-02-04 16:08

manager   ~0011840

No this does not demonstrate enthalpy conservation, it just shows the results are different, it doesn't even show which might be more "accurate". To test enthalpy conservation you need to integrate the enthalpy inlet and the outlet and compare the difference with the source.


2021-02-04 16:26

reporter   ~0011841

I totally agree with you, the temperature is not a conserved scalar.
I have checked the enthalpy balance and it is respected, but the calculated temperature is not in agreement with the expected one.
To me, there could be a problem in the conversion from enthalpy to temperature.


2021-02-04 17:01

manager   ~0011842

> To me, there could be a problem in the conversion from enthalpy to temperature.

Can you be more specific? What problem do you see in the conversion from enthalpy to temperature?


2021-02-04 18:04

manager   ~0011846

> in a more complex simulation, I have found a difference of 30°C

I don't have such a case to study.


2021-02-04 18:05

manager   ~0011847

> Also, it is not clear to me what do you mean with integrate the enthalpy.

You patch integrate the temperature flux, you should patch integrate the enthalpy flux.


2021-02-05 08:54

reporter   ~0011848

Dear Henry,

> You patch integrate the temperature flux, you should patch integrate the enthalpy flux.
thank you for your explanation, now it is clear

I have attached the outlet temperature profiles in what I have called "more complex case". This file shows the T comparison between the version OF6, OF7, OF8 taking into account the energy source.
> To test enthalpy conservation you need to integrate the enthalpy inlet and the outlet and compare the difference with the source.
I did, and in all cases, the difference of the integrated enthalpy between inlet and out is equal to the source. As you can see, the difference between OF and OF7/OF8 is around 30°C
I have also compared the outlet temperature between versions OF6 and OF8 without an energy source. In this case, the two versions are basically identical.


comparisons.pptx (707,584 bytes)


2021-02-05 09:25

manager   ~0011849

> I have attached the outlet temperature profiles in what I have called "more complex case".

Unfortunately this is not sufficient for me to be able to reproduce and analyse the difference.


2021-02-05 09:36

reporter   ~0011850

> Unfortunately this is not sufficient for me to be able to reproduce and analyse the difference.
 I know, but unfortunately, I cannot share this simulation, we don't have a non-disclosure agreement. For this reason, in the beginning, I have shared a simple case, which shows a small difference between the OF6 and OF8


2021-02-05 09:45

manager   ~0011851

Can you create a case you can share which shows the 30deg difference?


2021-02-05 09:57

reporter   ~0011852

Yes, I can.
It requires a little time to create a "similar" mesh and run the cases with OF6 and OF8. Moreover, I won't be able to share the results with this procedure due to the size of the case


2021-02-05 10:06

manager   ~0011853

If the mesh is generated with blockMesh then the case without the generated mesh would be small and you could share it. Without a test-case I would not be able to diagnose any problems.


2021-02-05 10:13

reporter   ~0011854

I am not able to generate this mesh via blockMesh. I have to create the mesh via Ansys and then to convert it in openfoam format via fluentMeshToFoam utility.


2021-02-05 10:17

manager   ~0011855

> I am not able to generate this mesh via blockMesh

Why not? blockMesh is supplied with OpenFOAM free of charge and open-source, you can use it.


2021-02-05 10:33

manager   ~0011857

If you wish you can contact CFD Direct to arrange a support contract and we can work with you directly.


2021-02-05 11:21

reporter   ~0011859

Dear Herny,
let me try to create a similar mesh via blockMesh


2021-02-05 11:32

manager   ~0011860

Does the mesh have to be similar to reproduce the problem? Does it have to be 3D? Does it have to have 7M cells?


2021-02-05 11:40

reporter   ~0011861

In order to have the same order of magnitude (30°C), I don't know if the mesh has to be in 3D, 7M of cells and so on.
The simple case which I sent yesterday reproduces this trend, even if the difference is smaller (3°C)


2021-02-05 12:23

manager   ~0011862

The 3°C difference is due to a change is the handling of the departure function for the incompressible temperature-dependent density.
It is not at all clear from what you have provided that the 30°C difference relates to this and without a simple test-case that reproduces THIS issue I cannot investigate further.


2021-02-05 12:28

reporter   ~0011863

I understand. Please, let me a few time to setup a new case via blockMesh to reproduce this behavior


2021-02-10 10:58

reporter   ~0011872

Dear Henry,
I have performed another simulation changing the thermo dictionary. Instead to apply heRhoThermo with polynomial density as done in OpenFoam version 6, I have applied hePsiThermo with PengRobinson eos with OpenFoam version 8. This simulation works and I obtained the expected results. Now the difference with OF6 is around 2°C.

I want to thank you for your patience and I really apologize for wasting your time.


Issue History

Date Modified Username Field Change
2021-02-03 14:19 tiziano.maffei New Issue
2021-02-03 15:00 wyldckat Note Added: 0011831
2021-02-03 15:06 tiziano.maffei File Added: TjunctionEnthalpySourceOF8.7z
2021-02-03 15:06 tiziano.maffei Note Added: 0011832
2021-02-04 14:17 henry Note Added: 0011833
2021-02-04 14:55 tiziano.maffei Note Added: 0011834
2021-02-04 15:04 henry Note Added: 0011835
2021-02-04 15:41 tiziano.maffei Note Added: 0011836
2021-02-04 15:48 henry Note Added: 0011837
2021-02-04 15:54 henry Note Added: 0011838
2021-02-04 16:02 tiziano.maffei Note Added: 0011839
2021-02-04 16:08 henry Note Added: 0011840
2021-02-04 16:26 tiziano.maffei Note Added: 0011841
2021-02-04 17:01 henry Note Added: 0011842
2021-02-04 18:04 henry Note Added: 0011846
2021-02-04 18:05 henry Note Added: 0011847
2021-02-05 08:54 tiziano.maffei File Added: comparisons.pptx
2021-02-05 08:54 tiziano.maffei Note Added: 0011848
2021-02-05 09:25 henry Note Added: 0011849
2021-02-05 09:36 tiziano.maffei Note Added: 0011850
2021-02-05 09:45 henry Note Added: 0011851
2021-02-05 09:57 tiziano.maffei Note Added: 0011852
2021-02-05 10:06 henry Note Added: 0011853
2021-02-05 10:13 tiziano.maffei Note Added: 0011854
2021-02-05 10:17 henry Note Added: 0011855
2021-02-05 10:33 henry Note Added: 0011857
2021-02-05 11:21 tiziano.maffei Note Added: 0011859
2021-02-05 11:32 henry Note Added: 0011860
2021-02-05 11:40 tiziano.maffei Note Added: 0011861
2021-02-05 12:23 henry Note Added: 0011862
2021-02-05 12:28 tiziano.maffei Note Added: 0011863
2021-02-10 10:58 tiziano.maffei Note Added: 0011872
2021-02-10 12:15 henry Assigned To => henry
2021-02-10 12:15 henry Status new => closed
2021-02-10 12:15 henry Resolution open => no change required