View Issue Details

IDProjectCategoryView StatusLast Update
0003046OpenFOAMBugpublic2018-08-21 18:58
ReporterjanGaertner Assigned Tohenry  
PrioritynormalSeverityminorReproducibilitysometimes
Status closedResolutionsuspended 
PlatformGNU/LinuxOSUbuntu OS Version14.04 LTS
Summary0003046: Initialization error of dgdt in createFields.H in compressibleInterFoam
DescriptionIn createFields.H the dgdt field in compressibleInterFoam is initialized with
dgdt = alpha1*fvc::div(phi)

However, for the first iteratrion the source term should be zero in the alpha equation:
(Sp*alpha1 + Su + alpha1*fvc::div(phi))
With,
dgdt = Sp*alpha1 + Su

This is the case when dgdt is initialized with:
dgdt = -alpha1*fvc::div(phi);

If the field is initialized incorrectly it leads to an error and to a mass imbalance of the solver.
Steps To ReproduceRun a case with a nozzle flow and a backfacing step and set the pressure as a fixed value boundary condition at inlet and outlet.
The mass will not be conserved and errors are in the range above 5%.
TagscompressibleInterFoam

Activities

henry

2018-08-16 15:35

manager   ~0009961

> However, for the first iteratrion the source term should be zero in the alpha equation:

Why? If div(phi) is non-zero then dgdt should be non-zero. Or do you think that dgdt should always be zero, only zero at the beginning of each time-step or zero only at the start of the run (time 0)?

henry

2018-08-16 15:49

manager   ~0009962

Is the issue that the initial phi field you have in your case is not conservative?
What do you get if you initialize your case with potentialFoam?

janGaertner

2018-08-17 08:34

reporter   ~0009966

The case is that if we start with a field of zero divergence, div(phi)=0. The terms on the RHS of the alpha equation should equal to zero as well.
If dgdt is initalized with dgdt = alpha1*fvc::div(phi) then this adds together with the divU field and the Su expression in the alpha equation. As a result the RHS is never zero, even for a divergence free field.

janGaertner

2018-08-17 08:36

reporter   ~0009967

Take a very simple test and run with one iteration and one outer loop with the compressibleInterFoam solver and see that the alpha field changes. Even if the initial divergence field is zero.

henry

2018-08-17 08:52

manager   ~0009968

Last edited: 2018-08-17 08:53

If div(phi)=0 then surely dgdt = alpha1*fvc::div(phi) is also 0?

It is not clear what change you want, can you clarify?

Do you think that dgdt should always be zero, only zero at the beginning of each time-step or zero only at the start of the run (time 0)?

Is the issue that the initial phi field you have in your case is not conservative?
What do you get if you initialize your case with potentialFoam?

janGaertner

2018-08-17 09:33

reporter   ~0009969

Ah yes, of course it is a bad example.

Case:
div(phi) not equal 0 in intial condition

AlphaEqn: (as implemented in VoF/alphaEqn.H)

ddt(alpha) + div(alpha*phi) = dgdt + alpha*div(phi)

Maybe I have an error how I think about this. Should it not be, that the RHS equals to zero in the first iteration and alpha only changes due to the divergence field on the LHS?

Currently it is implemented as:

ddt(alpha) + div(alpha*phi) = alpha*div(phi) + alpha*div(phi)

=> ddt(alpha) + alpha*div(phi) + phi*grad(alpha) = 2*alpha*div(phi)
=> ddt(alpha) + phi*grad(alpha) = alpha*div(phi)
=> DDt(alpha) = alpha*div(phi)

janGaertner

2018-08-17 09:52

reporter   ~0009970

One Note to add:
Assume a liquid vapor mix and liauid is represented with alpha1 and is nearly incompressible, so that rho1 approx. constant.

Then the alpha transport equation is:
ddt(rho1*alpha1) + div(rho1*alpha1*phi) = 0
or for nearly constant rho1:

=> ddt(alpha1) + div(alpha1*phi) = 0

This equation is only solved in alphaEqn.H if dgdt is initialized with: dgdt = -alpha1*div(phi)

henry

2018-08-17 10:19

manager   ~0009971

What do you get if you initialize your case with potentialFoam?

janGaertner

2018-08-17 10:35

reporter   ~0009972

not if I initialize with potential Foam.

But I will have cases with phase change and in these the divergence does not have to be zero.

Generally, I believe the descritpion as it is implemented is ok for most cases. But probably is not if the initial field is not well defined.

So changing the inital value is a small thing to guarantee that the equation system makes sense in these cases too.

henry

2018-08-17 11:03

manager   ~0009973

> not if I initialize with potential Foam.

That is what I get.

> But I will have cases with phase change and in these the divergence does not have to be zero.

OK, in which case dgdt should not be 0 either.

> Generally, I believe the descritpion as it is implemented is ok for most cases. But probably is not if the initial field is not well defined.

Right, and for these cases I would recommend using potentialFoam to make the initial field well defined.

> So changing the inital value is a small thing to guarantee that the equation system makes sense in these cases too.

What about restart?

henry

2018-08-17 11:46

manager   ~0009974

I have run some tests and confirmed that dgdt converges to alpha1*fvc::div(phi) as it should and this is correct initialization for dgdt. If your initial phi field is very non-conservative then there is a problem with this initialization for dgdt but this is not the only problem with the start-up. It would be best to ensure the initial phi field is consistent and appropriate, failing that you could use a large number of outer correctors at the beginning of the run to obtain convergence and conservation.

Have you tried increasing nOuterCorrectors? For my cases this resolves the problem.

janGaertner

2018-08-17 13:37

reporter   ~0009975

I do not understand why dgdt converges to alpha1*fvc::div(phi). Expecially if the phase 1 is nearly incompressible. As the RHS of the alpha equation is then always:
ddt(alpha1) + div(alpha*phi) = 2* alpha1 * fvc::div(phi)

henry

2018-08-21 18:58

manager   ~0009984

The proposed work-around will cause incorrect restart. It is not clear that any change is necessary and some important questions remain unanswered by the reporter. The general solution to the problem is to apply better initialization and/or more outer iterations early in the run to converge the initial fields.

Issue History

Date Modified Username Field Change
2018-08-16 14:41 janGaertner New Issue
2018-08-16 14:41 janGaertner Tag Attached: compressibleInterFoam
2018-08-16 15:35 henry Note Added: 0009961
2018-08-16 15:49 henry Note Added: 0009962
2018-08-17 08:34 janGaertner Note Added: 0009966
2018-08-17 08:36 janGaertner Note Added: 0009967
2018-08-17 08:52 henry Note Added: 0009968
2018-08-17 08:53 henry Note Edited: 0009968
2018-08-17 09:33 janGaertner Note Added: 0009969
2018-08-17 09:52 janGaertner Note Added: 0009970
2018-08-17 10:19 henry Note Added: 0009971
2018-08-17 10:35 janGaertner Note Added: 0009972
2018-08-17 11:03 henry Note Added: 0009973
2018-08-17 11:46 henry Note Added: 0009974
2018-08-17 13:37 janGaertner Note Added: 0009975
2018-08-21 18:58 henry Assigned To => henry
2018-08-21 18:58 henry Status new => closed
2018-08-21 18:58 henry Resolution open => suspended
2018-08-21 18:58 henry Note Added: 0009984