[Whyman]

Hi guys,

i have a problem with alpha1 in openFoam 1.6 ext. The value reduces with the time. Here is the output:



/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 1.6-ext
Exec : interGravitationalFoam
Date : Oct 16 2012
Time : 10:56:34
Host : Alderaan
PID : 3367
Case : /media/F2D48837D487FFD9/LAVORO/OPENFOAM/REXUS/test
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g

Loading Accelerations.data

Loading Omega.data
Calculating field g.h


Reading inertial Origin

Reading field pd

Reading field alpha1

Reading field U

Reading/calculating face flux field phi

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model kOmegaSST
kOmegaSSTCoeffs
{
alphaK1 0.85034;
alphaK2 1;
alphaOmega1 0.5;
alphaOmega2 0.85616;
gamma1 0.5532;
gamma2 0.4403;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
a1 0.31;
c1 10;
}

time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG: Solving for pcorr, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
Courant Number mean: 0 max: 0 velocity magnitude: 0
Omega1 initialized to (0 0 0)
Omega2 initialized to (0 0 0)

Starting time loop

Courant Number mean: 0 max: 0 velocity magnitude: 0
deltaT = 0.000595238
Time = 0.000595238

Calculating field g.h

Setting accelerations to (-112.644 1.09078 -0.109779)
Setting omega to (-0.00875158 -0.0815817 0.0820055)
Time at instant 1 is: 0
Time at instant 2 is: 0.000595238
Omega at instant 1 is: (0 0 0)
Omega at instant 2 is: (-0.00875158 -0.0815817 0.0820055)
Instantaneous angular acceleration is: (-14.7027 -137.057 137.769)
Previous Liquid phase volume fraction = 2.25683e-07 Min(alpha1) = 0 Max(alpha1) = 1
DILUPBiCG: Solving for alpha1: solution singularity
Liquid phase volume fraction = 2.25683e-07 Min(alpha1) = 0 Max(alpha1) = 1
DICPCG: Solving for pd, Initial residual = 1, Final residual = 9.62755e-08, No Iterations 329
time step continuity errors : sum local = 1.16224e-11, global = -1.58311e-23, cumulative = -1.58311e-23
time step continuity errors : sum local = 1.16224e-11, global = -1.58311e-23, cumulative = -3.16623e-23
DILUPBiCG: Solving for omega, Initial residual = 0.000655187, Final residual = 1.02048e-07, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 1.13822e-06, No Iterations 1
ExecutionTime = 27.63 s ClockTime = 27 s

Centre of Mass: (3.14552e-07 -4.44441e-05 0.000106692)

Courant Number mean: 0.000712776 max: 0.0760438 velocity magnitude: 10.4599
deltaT = 0.000705782
Time = 0.00130102

Calculating field g.h

Setting accelerations to (-113.489 1.10463 -0.123626)
Setting omega to (-0.00919504 -0.0653216 0.0650064)
Time at instant 1 is: 0.000595238
Time at instant 2 is: 0.00130102
Omega at instant 1 is: (-0.00875158 -0.0815817 0.0820055)
Omega at instant 2 is: (-0.00919504 -0.0653216 0.0650064)
Instantaneous angular acceleration is: (-0.628319 23.0383 -24.0855)
Previous Liquid phase volume fraction = 2.25683e-07 Min(alpha1) = 0 Max(alpha1) = 1
DILUPBiCG: Solving for alpha1, Initial residual = 0.00503263, Final residual = 5.05418e-09, No Iterations 2
Liquid phase volume fraction = 2.25683e-07 Min(alpha1) = 0 Max(alpha1) = 0.892477
DICPCG: Solving for pd, Initial residual = 0.918833, Final residual = 9.18764e-08, No Iterations 329
time step continuity errors : sum local = 1.97884e-11, global = 3.14334e-22, cumulative = 2.82672e-22
time step continuity errors : sum local = 1.97884e-11, global = 3.14334e-22, cumulative = 5.97006e-22
DILUPBiCG: Solving for omega, Initial residual = 0.000743617, Final residual = 2.82979e-08, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.333835, Final residual = 2.95525e-07, No Iterations 1
ExecutionTime = 46.73 s ClockTime = 47 s

Centre of Mass: (1.9406e-07 -3.56954e-05 8.55697e-05)

Courant Number mean: 0.000676682 max: 0.0735919 velocity magnitude: 8.53713
deltaT = 0.000839638
Time = 0.00214066

Calculating field g.h

Setting accelerations to (-114.494 1.1211 -0.1401)
Setting omega to (-0.0097226 -0.0459778 0.0447833)
Time at instant 1 is: 0.00130102
Time at instant 2 is: 0.00214066
Omega at instant 1 is: (-0.00919504 -0.0653216 0.0650064)
Omega at instant 2 is: (-0.0097226 -0.0459778 0.0447833)
Instantaneous angular acceleration is: (-0.628319 23.0383 -24.0855)
Previous Liquid phase volume fraction = 2.25683e-07 Min(alpha1) = 0 Max(alpha1) = 0.892477
DILUPBiCG: Solving for alpha1, Initial residual = 0.00309017, Final residual = 2.7438e-06, No Iterations 1
Liquid phase volume fraction = 2.25688e-07 Min(alpha1) = -1.11844e-97 Max(alpha1) = 0.824762
DICPCG: Solving for pd, Initial residual = 0.617329, Final residual = 9.87751e-08, No Iterations 287
time step continuity errors : sum local = 1.2306e-11, global = -1.7362e-24, cumulative = 5.95269e-22
time step continuity errors : sum local = 1.2306e-11, global = -1.7362e-24, cumulative = 5.93533e-22
DILUPBiCG: Solving for omega, Initial residual = 0.000872859, Final residual = 9.36405e-09, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.263603, Final residual = 1.13085e-07, No Iterations 1
ExecutionTime = 64.18 s ClockTime = 64 s

Centre of Mass: (1.11854e-07 -3.07893e-05 7.3643e-05)

Courant Number mean: 0.00056678 max: 0.0512143 velocity magnitude: 5.03397
deltaT = 0.00099707
Time = 0.00313773

Calculating field g.h

Setting accelerations to (-115.687 1.14066 -0.159662)
Setting omega to (-0.0103491 -0.0230069 0.0207683)
Time at instant 1 is: 0.00214066
Time at instant 2 is: 0.00313773
Omega at instant 1 is: (-0.0097226 -0.0459778 0.0447833)
Omega at instant 2 is: (-0.0103491 -0.0230069 0.0207683)
Instantaneous angular acceleration is: (-0.628319 23.0383 -24.0855)
Previous Liquid phase volume fraction = 2.25688e-07 Min(alpha1) = -1.11844e-97 Max(alpha1) = 0.824762
DILUPBiCG: Solving for alpha1, Initial residual = 0.00133358, Final residual = 3.93363e-07, No Iterations 1
Liquid phase volume fraction = 2.25688e-07 Min(alpha1) = -5.272e-08 Max(alpha1) = 0.79509



I have set buoyantPressure for p and pd and I run a modified interFoam.

Why do I have this result???

Please help me!!!

[akidess]

You don't give a lot of info, but try using explicit MULES for alpha1. Use a bounded divergence scheme and first order Euler.

[Whyman]

Which info do you need in particular?

I set MULES as Explicit, and Euler as time scheme. I don't know exactly what you mean with "bounded scheme"

Anyway it doesn't change, but still gives back the same error.

[akidess]

Helpful information would be the content of fvSchemes and fvSolution, specifying if you have any special boundary conditions (contact angles, inlets/outlets), or a problematic mesh.

About the bounded schemes - see "4.4.1.1 Schemes for strictly bounded scalar fields" in the user guide. Examples of such schemes are vanLeer01 and Gamma01.

- Anton

Originally Posted by Whyman

Which info do you need in particular?

I set MULES as Explicit, and Euler as time scheme. I don't know exactly what you mean with "bounded scheme"

Anyway it doesn't change, but still gives back the same error.

[Whyman]

Ok, I attach all the files for boundary condition plus fvSolution and fvScheme.

There is the contact angle acting in the domain.

[Tobi]

Originally Posted by Whyman

Ok, I attach all the files for boundary condition plus fvSolution and fvScheme.

There is the contact angle acting in the domain.

Hi Whyman (nice name by the way)



Your Problem description
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

I can not see any problem in your case. You said that the value of alpha1
goes below 1. Thats does not mean to be an error or an numerical problem. I had a look
in your time folder 0 and saw that there is no inlet for alpha1 that is
set to fixedValue.


I have got a notation that ...
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
you are init your alpha1 with the setFields command in a certain
region (like the damnBreak tutorial). While solving your case
the fraction 1 is distributed in the region (maybe due to gravitation)
and your fraction value reaches values below 1 due to distribution and the
fact that your fraction 1 is a thin layer or something like that.

Further more: You set your BC for U to fixedValue (0 0 0).
So I expect that you are in a closed box or sth. like that.


Experience
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
A year ago I made a simulation of a chocolate donot. Therefor I set a chocolate
field like a plate into my mesh and let it fall down on my donut due to gravitation.
I realized that the fraction is distributed due to the gravitation and velocity becouse
my chocolate plate was very thin and I used the visco of water and not of chocolate,
The results was, that my maxAlpha1 value goes down below 1.

That was logical couse the cells which belong to alpha1 = 1 take a high of 4 cells.
After the first few timesteps I got values below 1. I had a look at the first calculation
steps and realized that problem.


To do
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
First of all, write out the first time steps and have a look at it with paraview.
It should be possible to see your problem, or if your calculation is correct or wrong.

If everything is okay check your boundary conditions, initialization and so on. After
that you can go on and check your schemes. Like akidess told you use first order euler schemes
for time derivation and bounded scheme for alpha1.

But if the problem occure with MULES and DILUPBiCG it seems like the problem I mentioned and that
its not a problem in the numerical schemes.

I have no idea about your geometry and initialization but I think its due to distribution of the fraction
and the value goes under 1 couse your fraction 1 is a thin layer or sth. like that.


Further more
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
its possible that I skipped some details?.?.?
Without a concrete case I can not help you.
I grope on the dark.


But I Hope that this will help you a bit


Regard
Tobi

[akidess]

I agree with Tobi. I understood your original problem as that the volume integral of alpha1 is reducing (i.e. you are losing mass), but now I realized you meant that your maximum value is reducing (i.e. your field is dispersing). Follow the steps suggested by Tobi.

[Tobi]

Hi Anton,

i think your word "dispersed" is much common then my word used "distribution"

Tobi

[Whyman]

Hi guys,

thanks for your help....

Anyway, the domain is a cylindrical container, with inside some liquid and gas, subjected to acceleration varying with time.

Because I have modified a bit the interFoam solver, i send to you my UEqn and pEqn files: maybe i have made some mistakes in writing the code.

pEqn
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);

U = rUA*UEqn.H();

surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);

adjustPhi(phiU, U, pd);

phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();



for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rUAf, pd) == fvc::div(phi)
);

pdEqn.setReference(pdRefCell, pdRefValue);

if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
}

if (nonOrth == nNonOrthCorr)
{
phi -= pdEqn.flux();
}
}



U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();

# include "continuityErrs.H"

p = pd + rho*gh;

if (pd.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pdRefCell)
);
}

}


UEqn

surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
+ Fcentrifugal
+ Feuler
+ Fcoriolis
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf()))
);

UEqn.relax();

if (momentumPredictor)
{

solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}



I have also written the following lines in the main code, to take the gravity variation in account:

dimensionedVector acc("Acceleration", dimAcceleration,it(runTime.timeOutputValue()));
g = acc;

Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());


Maybe now my problem is clearer (I hope).