[rahulksoni]

I have found a typical problem with OpenFoam multiphase simulations. When trying to simulate with multiphaseInterFOAM, it runs well. But while doing multiphaseEulerFoam, it always results in shooting in Courant number and eventual failure.
I have 4 phases air, water, solid_phase1, solid_phase2. The last two phases having density 2600 and 4000 (higher than water). I have tried almost all possible combinations to get the simulation run but no luck. At the end (and following some reference over the internet), I can conclude that it is OpenFoam instability issue. Since there is a huge difference between the density of air and other phases. So, when momentum exchange occurs there is a velocity jump for air due to its very low density. It is noticeable that only momentum exchange calculation between phases at interface is extra in multiphaseEulerFoam when compared to multiphaseInterFoam

[BlnPhoenix]

Hello,

1) I'm not sure modelling solids with standard multiphaseeulerfoam is a valid approach. How do you set the viscosities?

2) From your post it is not clear if you're pointing to some form of bug or experiencing normal running problems. Your post is very unspecific regarding the setup and your problem while running. Is the solver crashing or is it just taking it's time to solve? One can only wonder...

[rahulksoni]

Originally Posted by BlnPhoenix

Hello,

1) I'm not sure modelling solids with standard multiphaseeulerfoam is a valid approach. How do you set the viscosities?

Estimating viscosity is difficult in such cases. Tried different models, but no luck. For time being (running only), doing it with liquid viscosity only.

2) From your post it is not clear if you're pointing to some form of bug or experiencing normal running problems. Your post is very unspecific regarding the setup and your problem while running. Is the solver crashing or is it just taking it's time to solve? One can only wonder...

I am not pointing to bug, it is a running problem, as far as I can understood from the error message. The Courant number shoots out. Following is the error message:

Courant Number mean: 0.000339307 max: 6375.3
deltaT = 1.45842e-17
Time = 0.01772

PIMPLE: iteration 1
MULES: Solving for alpha.air
air volume fraction, min, max = 0.999865 0 1.00031
MULES: Solving for alpha.water
water volume fraction, min, max = 0.000124458 0 0.92345
MULES: Solving for alpha.phase_light
phase_light volume fraction, min, max = 4.05842e-06 0 0.0318035
MULES: Solving for alpha.phase_dense
phase_dense volume fraction, min, max = 6.76404e-06 0 0.0511821
Phase-sum volume fraction, min, max = 1 1 1.00032
MULES: Solving for alpha.air
air volume fraction, min, max = 0.999865 0 1.00031
MULES: Solving for alpha.water
water volume fraction, min, max = 0.000124458 0 0.92345
MULES: Solving for alpha.phase_light
phase_light volume fraction, min, max = 4.05842e-06 0 0.0318035
MULES: Solving for alpha.phase_dense
phase_dense volume fraction, min, max = 6.76404e-06 0 0.0511821
Phase-sum volume fraction, min, max = 1 1 1.00032
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#5 Foam::dragModels::SyamlalOBrien::K(Foam::Geometric Field<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#6 Foam::multiphaseSystem::dragCoeffs() const at ??:?
#7 ? at ??:?
#8 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9 ? at ??:?

[BlnPhoenix]

Can you post CheckMesh output.

[rahulksoni]

Originally Posted by BlnPhoenix

Can you post CheckMesh output.

rial_from_simFlow2$ checkMesh
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 3.0.x-ac3f6c67e02f
Exec   : checkMesh
Date   : Oct 15 2017
Time   : 20:01:25
Host   : "rahul-HP-Z600-Workstation"
PID    : 14329
Case   : /media/rahul/Passport/OpenFOAM_Simulations/spiral/Attempt04/eulerFoam/mef_1.2L_Iron_trial_from_simFlow2
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

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

--> FOAM Warning : 
    From function dlOpen(const fileName&, const bool)
    in file POSIX.C at line 1179
    dlopen error : libturbulenceModelSchemes.so: cannot open shared object file: No such file or directory
--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName&, const bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 99
    could not load "libturbulenceModelSchemes.so"
Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           1312472
    faces:            2824487
    internal faces:   2262400
    cells:            764762
    faces per cell:   6.65159
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     612300
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     152462
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            6   41017
            7   34
            8   1664
            9   70808
           10   36
           11   927
           12   22454
           13   2
           14   4
           15   14406
           18   1108
           21   2

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points                  Surface topology
                  outlet     3713     5011  ok (non-closed singly connected)
                   inlet     1433     2039  ok (non-closed singly connected)
        atmospheric_face    86207   103207  ok (non-closed singly connected)
    spiralExceptAtmosphericPatch   470734   486954  ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.295262 -0.295262 0.5) (0.295262 0.295262 2.3)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (-1.94719e-17 -3.23215e-16 -1.81718e-14) OK.
    Max cell openness = 3.71735e-16 OK.
    Max aspect ratio = 4.65624 OK.
    Minimum face area = 1.10135e-07. Maximum face area = 0.00120969.  Face area magnitudes OK.
    Min volume = 1.13664e-10. Max volume = 2.7218e-05.  Total volume = 0.112692.  Cell volumes OK.
    Mesh non-orthogonality Max: 73.4373 average: 18.577
   *Number of severely non-orthogonal (> 70 degrees) faces: 939.
    Non-orthogonality check OK.
  <<Writing 939 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 3.93696 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

[BlnPhoenix]

So, altough it says "ok" in CheckMesh you should look at your severly non-orthogonal cells. MultiphaseEulerFoam doesn't like this and can crash. Also, try to find the Cells with high Courant number. If this doesn't help i would re-check Boundary Conditions.

[rahulksoni]

Originally Posted by BlnPhoenix

So, altough it says "ok" in CheckMesh you should look at your severly non-orthogonal cells. MultiphaseEulerFoam doesn't like this and can crash. Also, try to find the Cells with high Courant number. If this doesn't help i would re-check Boundary Conditions.

Thanks a lot BlnPhoenix for the hints. I was guessing to check the individual meshes for their orthogonality and for high-courant number. But How would I do that in paraFoam (or if any other option is there) because when I try to clip or slice the system, it takes a lot more time as number of meshes are too high. It takes a long time for each refresh. Is there any way to filter cells in paraFoam to see only those cells with high non-orthogonality or courant number/velocity. I couldn't find such an option or may be difficult to implement.

The boundary condition folder for your convenience is uploaded.

[BlnPhoenix]

When you used the "Co" Tool, can you try to use a Contour Plot for Co to locate the Cells?


Edit: Since you asked about the Non-Ortho Cells also. I know that after CheckMesh there should be a folder with the Cells that have high non-ortho value. There should be a way to visualize them, can not recall it exactly ho it's done atm though. Try to search the forum for .vtk and visualization. I'm sure its on here somewhere.

I can not check your .tar.gz file because i'm in windows.

[rahulksoni]

Originally Posted by BlnPhoenix

When you used the "Co" Tool, can you try to use a Contour Plot for Co to locate the Cells?


Edit: Since you asked about the Non-Ortho Cells also. I know that after CheckMesh there should be a folder with the Cells that have high non-ortho value. There should be a way to visualize them, can not recall it exactly ho it's done atm though. Try to search the forum for .vtk and visualization. I'm sure its on here somewhere.

I can not check your .tar.gz file because i'm in windows.

Thanks Thanks for the wonderful information on "Co" utility that I was totally unaware. Combined with the help from somewhere else on paraview Threshold filter (How to choose the right courant number? I can produce the attached snapshot of my case.

Certainly, there are regions of high velocities (top-right image) and regions of high courant number (bottom-right image). But there is a mismatch in region selection. Not all regions/cells of high velocity possess the high courant number, indicating corresponding bigger cells. Thus, Courant number would be the true indicator for the problematic cells.
The second image shows the zoomed view of the problematic region. As obvious there are cells with non-orthogonality. But the question is how to get rid of them.

After analyzing the snapshots following things come to my mind:
1. Is problem really because of non-orthogonal cells. Saying the same because the cells with high courant number are all hex (as list from the obvious visibility).
2. Are courant number high because these cells are smaller than their neighbors?


PS: I have updated a .zip format of 0 folder in the previous post.

[rahulksoni]

Originally Posted by rahulksoni

Thanks Thanks for the wonderful information on "Co" utility that I was totally unaware. Combined with the help from somewhere else on paraview Threshold filter (How to choose the right courant number? I can produce the attached snapshot of my case.

Certainly, there are regions of high velocities (top-right image) and regions of high courant number (bottom-right image). But there is a mismatch in region selection. Not all regions/cells of high velocity possess the high courant number, indicating corresponding bigger cells. Thus, Courant number would be the true indicator for the problematic cells.
The second image shows the zoomed view of the problematic region. As obvious there are cells with non-orthogonality. But the question is how to get rid of them.

After analyzing the snapshots following things come to my mind:
1. Is problem really because of non-orthogonal cells. Saying the same because the cells with high courant number are all hex (as list from the obvious visibility).
2. Are courant number high because these cells are smaller than their neighbors?
3. The bigger question is why the flow as high velocity near the edges only? The top cover of the whole spiral is a uniform atmospheric patch. Then why should flow is high near the edges only? The air can outlet from the whole patch. Similarly, the whole circular face at the top-most surface is defined as inlet with proper inlet condition then how come again the velocities are high just near the edges. This is surprising for me.

PS: I have updated a .zip format of 0 folder in the previous post.

Am I missing something or my boundary conditions are wrong (as I am guessing from the flow pattern)? May be pressure is high in the main flow spiral region and that is why things are short-circuiting from inlet to atmospheric patch. Am I getting the wrong interpretation of term p_rgh????

The essence ( the Complete list can be found in attached 0.zip folder) of boundary conditions of my case is presented below:

p_rgh

internalField	uniform 0.0;

boundaryField
{
	inlet
	{
		type	fixedFluxPressure;
		value 	uniform 0.0;
	}
	outlet
	{
		type    fixedValue; 
                value   uniform 0;
	}
	atmospheric_face
	{
		type    fixedValue;
               value   uniform 0;
	}	
	spiralExceptAtmosphericPatch
	{
		type	zeroGradient; 
	}
}

phase velocity

internalField	uniform (0.0 0.0 0.0);

boundaryField
{
	outlet
	{
		type	pressureInletOutletVelocity;
		value	uniform (0.0 0.0 0.0);
	}
	atmospheric_face
	{
		type	pressureInletOutletVelocity;
		value	uniform (0.0 0.0 0.0);
	}
	inlet
	{
		type	fixedValue;
		value	uniform (-0.04217 -0.05622 -0.05622);
	}
	spiralExceptAtmosphericPatch
	{
		type	fixedValue;
		value	uniform (0.0 0.0 0.0);
	}
}

phase fraction

internalField	uniform 0.0;

boundaryField
{
	outlet
	{
		type		inletOutlet;
		inletValue	uniform 0.0;
	}
	atmospheric_face
	{
		type		inletOutlet;
		inletValue	uniform 0.0;
	}
	inlet
	{
		type		fixedValue;
		value		uniform 0.92;
	}
	spiralExceptAtmosphericPatch
	{
		type		zeroGradient;
	}
}

[BlnPhoenix]

From your BC's i can not make out anything obviously wrong.
So my guess is still the mesh. You could try non-orthogonal corrections. If you tried that already or it's not working, there is usually only one way - re-meshing. What is your meshing tool?

[BlnPhoenix]

Maybe you can also try to make the athmosphere as a wall not as an outlet. Just as a test, if this runs well.

[rahulksoni]

Originally Posted by BlnPhoenix

Maybe you can also try to make the athmosphere as a wall not as an outlet. Just as a test, if this runs well.

Tried but no luck.

[rahulksoni]

Originally Posted by BlnPhoenix

From your BC's i can not make out anything obviously wrong.
So my guess is still the mesh. You could try non-orthogonal corrections. If you tried that already or it's not working, there is usually only one way - re-meshing. What is your meshing tool?

I believe increasing the nNonOrthogonalCorrectors values in fvSolutions file is going to handle the orthogonal correcrions.

Second thing is that I have closely looked into the meshing in the region where courant number is high. I can see that the meshes are really fine there (as visible in the attached snapshot (left image)). It is because of the snapping control feature that tries to smooth the joints or edges by breaking down in increased number of cells. As top surface is atmospheric patch (just to give natural condition and for air outlet, if it is there), it is not of much concern and coarser mesh near that would be completely fine. But how do I achieve it. Bottom I am pasting my snappyHexMeshDict content. What condition is govering it.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
	version	2.0;
	class	dictionary;
	format	ascii;
	location	"system";
	object	snappyHexMeshDict;
}
castellatedMesh			true;
snap					true;
addLayers				false;
geometry
{
	atmospheric_face.stl
	{
		type			triSurfaceMesh;
		simflowType		stl;
		name			atmospheric_face;
		includedAngle	120.0;
	}
	spiralExceptAtmosphericPatch.stl
	{
		type			triSurfaceMesh;
		simflowType		stl;
		name			spiralExceptAtmosphericPatch;
		includedAngle	120.0;
	}
}
castellatedMeshControls
{
	locationInMesh		(0.05 -0.05 2.2);
	
	refinementSurfaces
	{	
		atmospheric_face
		{
			level		(1 1);
		}
		spiralExceptAtmosphericPatch
		{
			level		(2 4);				// These are min max values, refinement can go up to max if surface curvature is more than angle defined by keyword resolveFeatureAngle. See page 32 of https://openfoamwiki.net/images/f/f0/Final-AndrewJacksonSlidesOFW7.pdf
		}
	}
	refinementRegions
    {
        spiralExceptAtmosphericPatch
        {
            mode        distance;
            levels      ((0.010 3) (0.030 2) (0.050 1));//((0.005 3) (0.015 2) (0.03 1));122556 nCells
        }
        atmospheric_face
        {
            mode        distance;
            levels      ((0 0));
        }
    }
	limitRegions
	{
	}
	features				// Basically important for capuring the edges properly, not having much effect in current case
	(
		{
    		file "atmospheric_face.eMesh";
    		level 1;
  		}
  		{
    		file "spiralExceptAtmosphericPatch.eMesh";
    		level 5;
  		}
	);
	maxGlobalCells				1000000;
	maxLocalCells				100000;
	resolveFeatureAngle			30.0;//100.0;
	maxLoadUnbalance			0.1;
	nCellsBetweenLevels			1;
	gapLevelIncrement			0;
	minRefinementCells			3;
	allowFreeStandingZoneFaces	false;
}

meshQualityControls
{
	maxConcave					120.0;
	maxBoundarySkewness			5.0;//10.0;
	minFaceWeight				0.05;
	maxInternalSkewness			5.0;//9.0;
	minTwist					0.05;
	maxNonOrtho					0.0;//85.0;
	minDeterminant				0.01;
	minVolRatio					0.01;
	nSmoothScale				10;
	minTriangleTwist			-1.0;
	minVol						1.0E-20;
	minTetQuality				1.0E-20;
	minArea						-1.0;
	errorReduction				0.9;
	relaxed
	{
		maxNonOrtho				5.0;///75.0;
	}
}
mergeTolerance					1.0E-06;//1.0E-4;
debug							0;
combinePatches					true;

snapControls
{
	implicitFeatureSnap			false;
	nSmoothPatch				5;
	multiRegionFeatureSnap		false;
	nFeatureSnapIter			5;
	nSolveIter					50;
	explicitFeatureSnap			false;
	tolerance					0.5;
	nRelaxIter					5;
}

addLayersControls
{
	layers
	{
		        atmospheric_face
                {
                        nSurfaceLayers  		1;
                        expansionRatio  		1.25;
                        firstLayerThickness     0.2;
                }
                atmospheric_face_patch0
                {
                        nSurfaceLayers  		1;
                        expansionRatio  		1.25;
                        firstLayerThickness     0.2;
                }
                spiralExceptAtmosphericPatch
                {
                        nSurfaceLayers  		2;
                        expansionRatio  		1.25;
                        firstLayerThickness     0.2;
                }
                spiralExceptAtmosphericPatch_patch0
                {
                        nSurfaceLayers  		2;
                        expansionRatio  		1.25;
                        firstLayerThickness     0.2;
                }
	}
	nMedialAxisIter				100;
	minThickness				0.05;
	nSmoothNormals				3;
	maxThicknessToMedialRatio	0.5;
	nSmoothThickness			10;
	nSmoothSurfaceNormals		5;
	slipFeatureAngle			30.0;
	minMedianAxisAngle			90.0;
	nRelaxIter					5;
	concaveAngle				90.0;
	nGrow						0;
	nBufferCellsNoExtrude		0;
	nRelaxedIter				20;
	featureAngle				45.0;
	maxFaceThicknessRatio		0.5;
	nLayerIter					50;
	relativeSizes				true;
	firstLayerThickness			0.5;
	expansionRatio				1.25;
	layerTerminationAngle		-180.0;
	meshShrinker				displacementMedialAxis;
}

[BlnPhoenix]

Originally Posted by rahulksoni

Tried but no luck.

Ok.

Search for nNonOrthogonalCorrectors. You can set this value in fvsolution. You can try to run it with a value 2-3 if i'm not wrong. But this usually doesn't help in situations where non-orthogonality is very high.

I'm no expert in snappyhexmesh but from the glance of the dict i would try a value of 45 for maxNonortho. A value of 0 or 5 is not doable for snappy! Also a value of 85 is to high.
Also please try to set either implicitFeatureSnap or explicitFeatureSnap to "true" or vice versa.

[rahulksoni]

Finally, after a lot of hit and trial and help on CO utility to produce courant number, I found that it is due to improper meshing only. Actually, snappyhexMesh tries to create the fine mesh at the corners and near the edges for their proper shape recognition. Fine meshes have smaller cell length (d). And, as we know the Courant number is:

Co = U x t / d

where, U is fluid velocity in the cell, t is the current timestep size and d as said is cell length. being in the denominator d creates a lot of problems. A smaller value shoots the value of Co, forcing OpenFoam to reduce the timestep size (t) for manageable converged solutions.

So, the only solution was to recreate the mesh with minimal fine meshes. Moreover, the regions with fine meshes that were creating problem were not of my interest of calculation. So, I thought to recreate the mesh.

As we can see from the attached first pic that Courant number is high on and near the side wall. And, it was because that side wall and bottom spiral wall were unibody causing the same snappy condition of fine mesh near the wall to apply for both of them, as seen in the left part of the second attached pic. Later, I have separated these two walls as different bodies and removed fine meshing condition near the side wall. The resulted improved mesh is shown on the right side of the second pic.

[Dyrney]

Dear,

I would like to know if someone could help me with how to set up the "phaseProperties" file (in "multiPhaseEulerFoam" solver) for more than one solid phase (solid_1-solid_2-air).
How to specify the continuous phase? How to set up the drag force and heat transfer?

Thanks.