[king1989]

Hi everyone,

I've been having trouble with OpenFOAM (buoyantPimpleFoam) for the past several weeks trying to simulate a tubular device with a rotating outer wall boundary condition. It looks something like this:
http://i.imgur.com/fCILHqU.png

There's an inlet at the bottom, an outlet at the top, and a cylinder in the middle. This cylinder is supposed to 'fake' an air-water boundary because I don't want to go through the trouble of running a multiphase simulation (at the moment). As you might have guessed, this device is a tubular centrifuge.

When my walls don't rotate, the velocity profile looks fine: maximum axial velocity occurs at the interface (where I have a free-slip boundary condition). However, when it does rotate, this happens:

http://i.imgur.com/g6S0xeL.gifv

The axial velocity starts out fine, but then there's these weird spots that develop where the velocity is unnaturally large and other spots where it's going in the opposite direction! What's even weirder is that when I look at the tangential components (ux and uy), they don't budge at all, so the rotation at the wall doesn't even propagate to the rest of the fluid. For reference, here are my U and p_rgh files.

P_rgh:

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 1e5;

boundaryField
{

inlet
{
    type            fixedFluxPressure;
}

outlet
{
    type            fixedValue;
    value           uniform 1e5;
}

walls
{
    type            fixedFluxPressure;
}

interface
{
type           zeroGradient;
}

}

U:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{

inlet
{
type        fixedValue;
value       uniform (0 0 0.1);
}   

outlet
{
type        pressureInletVelocity;
value       uniform (0 0 0);
}

walls
{
    type            rotatingWallVelocity;
    origin          (0 0 0);
    axis            (0 0 1);
    omega           39.79; //380 rpm
}

interface
{
type        slip;
}

}

I should also mention that changing to MRF doesn't help either. The solutions just straight up blow up and have big trouble converging. Given my inexperience with OpenFOAM, it's possible that the manner in which I implement MRF is completely wrong. Please let me know if you need any files from me.

[king1989]

Anyone....?

[blais.bruno]

From your picture, you seem to be using a tetrahedral mesh?

Have you put non-orthogonal correctors in your solver? I would suggest making a structured hexahedral mesh of your problem or at least a hex-dominant mesh, since tet mesh perform generally poorly in OpenFOAM.

Originally Posted by king1989

Anyone....?

[king1989]

Yes, my mesh is tetrahedral. When I check my fvSolution, it's something like this. So no, I'm not using any non-orthogonal correctors right now. Do you suggest using them? Sorry, I'm a bit of a noob at this and OpenFOAM's documentation is rather sparse and all over the place.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-9;
        relTol          0;
        smoother        DICGaussSeidel;
	nPreSweeps	0;
	nPostSweeps	2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
	maxIter		50;
    }

    pFinal
    {
        $p;
    }

    "U.*"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         3;
        tolerance       1e-09;
        relTol          1e-4;
    }

    "rho.*"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       0;
        relTol          0;
    }

    p_rgh
    {
        $p;
    }

    p_rghFinal
    {
        $p;
    }
    
    hFinal
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }
    h
    {
	$hFinal;
    }
}

PIMPLE
{
    transonic       no;
    nOuterCorrectors 100;
    nCorrectors     1;
    nNonOrthogonalCorrectors 0;
    momentumPredictor yes;
    pRefCell	0;
    pRefValue	0;

    residualControl
	{
		U
		{
			tolerance 1e-5;
			relTol	0;
		}
		p_rgh
		{
			tolerance 1e-4;
			relTol  0;
		}
	}
}

relaxationFactors
{
    fields
    {
	p_rgh	0.3;
	p_rghFinal  1;
    }
    equations
    {
	U	0.3;
	UFinal	1;
    }
}

// ************************************************************************* //

[king1989]

To add to the previous comment, the software I'm using (Salome) doesn't have an option for a hexahedral mesh. Is it imperative to switch to hexahedral or is there a workaround with a tetrahedral mesh?

[blais.bruno]

Salome can easily do hexahedral mesh, you just need to structure your mesh a bit. Obviously it won't be automatic, but generating a finite volume mesh is never really an automatic process. You can use snappyHexMesh of cfMesh also if you want, those are two good tools. In all cases, OpenFOAM performs poorly with tetrahedral meshes.

You need to put nNonOrthogonalCorrector to 2 or higher.

Cheers

Originally Posted by king1989

To add to the previous comment, the software I'm using (Salome) doesn't have an option for a hexahedral mesh. Is it imperative to switch to hexahedral or is there a workaround with a tetrahedral mesh?

[king1989]

Originally Posted by blais.bruno

Salome can easily do hexahedral mesh, you just need to structure your mesh a bit. Obviously it won't be automatic, but generating a finite volume mesh is never really an automatic process. You can use snappyHexMesh of cfMesh also if you want, those are two good tools. In all cases, OpenFOAM performs poorly with tetrahedral meshes.

You need to put nNonOrthogonalCorrector to 2 or higher.

Cheers

Thank you for your help. What do you mean by 'structure your mesh'? Are there specific rules I should follow when creating the geometry file in python for Salome to use (shown below for my current mesh):

# Creating Last Design

# Necessary for all geometries created by python code
import salome
salome.salome_init()
import GEOM
import math
from salome.geom import geomBuilder
geompy = geomBuilder.New(salome.myStudy)
gg = salome.ImportComponentGUI("GEOM")

# Creating -OZ vector
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
OZ_neg = geompy.MakeVectorDXDYDZ(0, 0, -1)

# -------------------- Parameters of the Geometry -------------------
# Points
# Cylinder
h_cyl = 0.2
r_cyl  = 0.03 

r_inlet = 0.005
r_outlet = 0.01
h_inletcyl = 0.001 
h_outletcyl = 0.001 

# Blades
z1_blade = 0.005
z2_blade = 0.195
blades_len = 0.029
blades_thick = 0.001




# ------------------- Creating Geometry ------------------------
bp = geompy.MakeVertex(0, 0, 0) # Base Point
bp_interface = geompy.MakeVertex(0, 0, 0.9 * z1_blade) # Base Point
bp_outlet = geompy.MakeVertex(0, 0, h_cyl) # Top point of the cylinder
bp_blades = geompy.MakeVertex(0, 0, z1_blade) # Top point of the cylinder


# Inlet, Outlet, and Main cylinder
cyl_inlet = geompy.MakeCylinder(bp, OZ_neg, r_inlet, h_inletcyl)  # Inlet cylinder
cyl_outlet = geompy.MakeCylinder(bp_outlet, OZ, r_outlet, h_outletcyl)  # Outlet cylinder
cyl = geompy.MakeCylinder(bp, OZ, r_cyl, h_cyl)  # Main cylinder

# Blades 
pt1 = geompy.MakeVertex(blades_thick, -blades_len, z1_blade)
pt2 = geompy.MakeVertex(-blades_thick, blades_len, z2_blade)
blade1 = geompy.MakeBoxTwoPnt(pt1, pt2)
angle = 90 * math.pi / 180
blade2 = geompy.MakeRotation(blade1, OZ, angle)

h_interface = h_cyl - 0.9*z1_blade
interface = geompy.MakeCylinder(bp_interface, OZ, 8.582815702e-3, h_interface)  # Air-water interface

# Partition and Cut
partition = geompy.MakePartition([cyl],[interface,cyl_inlet,cyl_outlet])
cut1 = geompy.MakeCut(partition,interface)

# -------------- Create Groups from cut1 ---------------------
inlet = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(inlet, [25])

outlet = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(outlet, [20])

walls = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(walls, [3,10,15])

interface = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(interface, [27,32])


 #------------ add objects in the study ----------------------

id_cut1 = geompy.addToStudy(cut1,"Cut1")
id_inlet = geompy.addToStudyInFather(cut1, inlet, "INLET")
id_outlet = geompy.addToStudyInFather(cut1, outlet, "OUTLET")
id_walls = geompy.addToStudyInFather(cut1, walls, "WALLS")
id_interface = geompy.addToStudyInFather(cut1, interface, "INTERFACE")

[blais.bruno]

What I mean by strucuture is to build a semi-structured mesh (combine different blocks, see where things can be structured). Salome does not have an automatic hexahedral mesher, but if you do the geometry block by block you can get an hex-dominant mesh.

Originally Posted by king1989

Thank you for your help. What do you mean by 'structure your mesh'? Are there specific rules I should follow when creating the geometry file in python for Salome to use (shown below for my current mesh):

# Creating Last Design

# Necessary for all geometries created by python code
import salome
salome.salome_init()
import GEOM
import math
from salome.geom import geomBuilder
geompy = geomBuilder.New(salome.myStudy)
gg = salome.ImportComponentGUI("GEOM")

# Creating -OZ vector
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
OZ_neg = geompy.MakeVectorDXDYDZ(0, 0, -1)

# -------------------- Parameters of the Geometry -------------------
# Points
# Cylinder
h_cyl = 0.2
r_cyl  = 0.03 

r_inlet = 0.005
r_outlet = 0.01
h_inletcyl = 0.001 
h_outletcyl = 0.001 

# Blades
z1_blade = 0.005
z2_blade = 0.195
blades_len = 0.029
blades_thick = 0.001




# ------------------- Creating Geometry ------------------------
bp = geompy.MakeVertex(0, 0, 0) # Base Point
bp_interface = geompy.MakeVertex(0, 0, 0.9 * z1_blade) # Base Point
bp_outlet = geompy.MakeVertex(0, 0, h_cyl) # Top point of the cylinder
bp_blades = geompy.MakeVertex(0, 0, z1_blade) # Top point of the cylinder


# Inlet, Outlet, and Main cylinder
cyl_inlet = geompy.MakeCylinder(bp, OZ_neg, r_inlet, h_inletcyl)  # Inlet cylinder
cyl_outlet = geompy.MakeCylinder(bp_outlet, OZ, r_outlet, h_outletcyl)  # Outlet cylinder
cyl = geompy.MakeCylinder(bp, OZ, r_cyl, h_cyl)  # Main cylinder

# Blades 
pt1 = geompy.MakeVertex(blades_thick, -blades_len, z1_blade)
pt2 = geompy.MakeVertex(-blades_thick, blades_len, z2_blade)
blade1 = geompy.MakeBoxTwoPnt(pt1, pt2)
angle = 90 * math.pi / 180
blade2 = geompy.MakeRotation(blade1, OZ, angle)

h_interface = h_cyl - 0.9*z1_blade
interface = geompy.MakeCylinder(bp_interface, OZ, 8.582815702e-3, h_interface)  # Air-water interface

# Partition and Cut
partition = geompy.MakePartition([cyl],[interface,cyl_inlet,cyl_outlet])
cut1 = geompy.MakeCut(partition,interface)

# -------------- Create Groups from cut1 ---------------------
inlet = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(inlet, [25])

outlet = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(outlet, [20])

walls = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(walls, [3,10,15])

interface = geompy.CreateGroup(cut1, geompy.ShapeType["FACE"])
geompy.UnionIDs(interface, [27,32])


 #------------ add objects in the study ----------------------

id_cut1 = geompy.addToStudy(cut1,"Cut1")
id_inlet = geompy.addToStudyInFather(cut1, inlet, "INLET")
id_outlet = geompy.addToStudyInFather(cut1, outlet, "OUTLET")
id_walls = geompy.addToStudyInFather(cut1, walls, "WALLS")
id_interface = geompy.addToStudyInFather(cut1, interface, "INTERFACE")