View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002073OpenFOAMBugpublic2016-04-29 10:322016-05-03 14:58
Reporterdemichie Assigned Tohenry  
PrioritynormalSeverityminorReproducibilitysometimes
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Fixed in Versiondev 
Summary0002073: density problem in reactingTwoPhaseEulerFoam
DescriptionIn several tests I have done with reactingTwoPhaseEulerFoam I have encountered a problem with negative values of the density of one of the phases, in particular with gas. It seems that the problem originates in pU/pEqn.H (or pUf/pEqn.H) here:

rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);

If I check the new value of rho1 and I compare it with psi1*p there is an increasing difference. This is due to the fact that (p_rgh - p_rgh_0) is not changing because of changes in pressure only, but also because rho is changing. It is updated with fluid.correctThermo() in EEqns.H.

If fluid.correctThermo() is used also in pEqn.H instead of the original lines

//rho1 += psi1*(p_rgh - p_rgh_0);
//rho2 += psi2*(p_rgh - p_rgh_0);

fluid.correctThermo()


the problem is solved and there are no more negative densities.
The same result is obtained introducing a field p0 before the pressure loop:

// Cache p prior to solve for density update
    volScalarField p_rgh_0("p_rgh_0", p_rgh);
    volScalarField p0("p0", p);


and then updating the densities with the change in pressure:


    //rho1 += psi1*(p_rgh - p_rgh_0);
    //rho2 += psi2*(p_rgh - p_rgh_0);
    rho1 += psi1*(p - p0);
    rho2 += psi2*(p - p0);





TagsNo tags attached.

Activities

henry

2016-04-29 10:53

manager   ~0006183

The problem with your approach is that phase continuity may be violated because the density update no longer corresponds to the p_rgh changes.

henry

2016-04-29 10:54

manager   ~0006184

Could you also test your cases in OpenFOAM-dev in which there have been many improvements to the multiphase solvers?

demichie

2016-04-29 13:33

reporter   ~0006185

I am compiling OpenFOAM-dev and I will test it later.

I have also tried another thing.
After rho is changed (in the main file with fluid.correct and in EEqns.H with fluid.correctThermo) now I am updating also p_rgh with:

p_rgh -= gh*(rho-rho_0);

where rho_0 is a value cached before the density update.

demichie

2016-04-29 14:01

reporter   ~0006186

I have just checked and the same problem appears with latest OpenFoam-dev.

henry

2016-04-29 14:22

manager   ~0006187

> p_rgh -= gh*(rho-rho_0);

will change p_rgh to be inconsistent with the p_rgh equation.

henry

2016-04-29 17:26

manager   ~0006196

Your proposals will violate phase continuity and not wise. Unfortunately I do not have any cases which reproduce the problem so cannot test any alternative approches to pressure and/or density clipping.

demichie

2016-04-29 17:40

reporter   ~0006197

Ok, but it is not clear to me why density should be updated with

rho1 += psi1*(p_rgh - p_rgh_0);


and not with

rho1 += psi1*(p - p0);

The last one is consistent with the equation of state, while using the second one the density of each phase changes also because of changes in the density of the mixture (affecting p_rgh).

In any case, I can try to create a simple case where the problem is present.
Unfortunately I will be away all the next week, so it is probable that I can upload the files the following week.

henry

2016-04-29 18:13

manager   ~0006198

rho1 += psi1*(p_rgh - p_rgh_0);

Is consistent with the linearization in the p_rgh equation which uses

correction(fvm::ddt(p_rgh));

to handle the variation of density due to pressure around the current state.

rho1 += psi1*(p - p0);

should be the same as

rho1 += psi1*(p_rgh - p_rgh_0);

provided p0 is stored at the same point p_rgh_0 is currently stored. Where in the code are you storing p0?

demichie

2016-04-29 18:29

reporter   ~0006199

I am storing p0 and p_rgh_0 at the same point:


    // Cache p prior to solve for density update
    volScalarField p_rgh_0("p_rgh_0", p_rgh);
    volScalarField p_0("p_0", p);

You can check that with the following two updates of the densities there are different results:

    // Update densities from change in p_rgh
    // rho1 += psi1*(p_rgh - p_rgh_0);
    // rho2 += psi2*(p_rgh - p_rgh_0);
    rho1 += psi1*(p - p_0);
    rho2 += psi2*(p - p_0);

henry

2016-04-29 18:39

manager   ~0006200

Interesting; if you can provide a suitable test case I will investigate further. Most likely the issue is to do with the statement

    p_rgh = p - rho*gh;

not updating the p_rgh boundary conditions. Do you find that the density unboundedness is on or near the boundary?

demichie

2016-04-30 18:47

reporter  

TestReacting.zip (1,012,668 bytes)

demichie

2016-04-30 18:50

reporter   ~0006203

I have uploaded a test case where the solver crash with this messages:

#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::sqrt(Foam::Field<double>&, Foam::UList<double> const&) at ??:?
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::sqrt<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#5 Foam::massTransferModels::Frossling::K() const at ??:?
#6 Foam::BlendedInterfacialModel<Foam::massTransferModel>::K() const at ??:?
#7 Foam::InterfaceCompositionPhaseChangePhaseSystem<Foam::MomentumTransferPhaseSystem<Foam::twoPhaseSystem> >::massTransfer() const at ??:?
#8 ? at ??:?
#9 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10 ? at ??:?
Floating point exception (core dumped)

As I wrote, the problem originates from a negative density for rho1.

henry

2016-04-30 18:55

manager   ~0006204

Do you find that the density unboundedness is on or near the boundary?

demichie

2016-04-30 18:57

reporter   ~0006205

I know that it happens also in the internal values, but I don't know if it is near the boundary.

henry

2016-04-30 19:25

manager   ~0006206

Have you tried thresholding the cells in which the density first goes negative to see if they are on the boundary?

demichie

2016-04-30 19:28

reporter   ~0006207

Unfortunately now I am on my laptop without OpenFOAM installed. I have to wait a couple of weeks when I will back in my office for additional checks.

henry

2016-05-03 14:57

manager   ~0006221

Thanks for providing an excellent test-case to demonstrate the issue. The fix you proposed only resolves part of the issue which relates to the need to update p_rgh after the density update at the beginning of pEqn.H:

commit 0f1921787e196ad772095b629760b18dc0d74b94

I will make the same change to OpenFOAM-3.0.x shortly.

Issue History

Date Modified Username Field Change
2016-04-29 10:32 demichie New Issue
2016-04-29 10:53 henry Note Added: 0006183
2016-04-29 10:54 henry Note Added: 0006184
2016-04-29 13:33 demichie Note Added: 0006185
2016-04-29 14:01 demichie Note Added: 0006186
2016-04-29 14:22 henry Note Added: 0006187
2016-04-29 17:26 henry Note Added: 0006196
2016-04-29 17:40 demichie Note Added: 0006197
2016-04-29 18:13 henry Note Added: 0006198
2016-04-29 18:29 demichie Note Added: 0006199
2016-04-29 18:39 henry Note Added: 0006200
2016-04-30 18:47 demichie File Added: TestReacting.zip
2016-04-30 18:50 demichie Note Added: 0006203
2016-04-30 18:55 henry Note Added: 0006204
2016-04-30 18:57 demichie Note Added: 0006205
2016-04-30 19:25 henry Note Added: 0006206
2016-04-30 19:28 demichie Note Added: 0006207
2016-05-03 14:57 henry Note Added: 0006221
2016-05-03 14:57 henry Status new => resolved
2016-05-03 14:57 henry Fixed in Version => dev
2016-05-03 14:57 henry Resolution open => fixed
2016-05-03 14:57 henry Assigned To => henry