segunda-feira, 10 de junho de 2013

How to set Icofoam to solve the Annular Flow problem

As you observed, I solved the problem already. This post is a step-by-step to help you to solve the same problem.

What do we need?

  1. Prepare a set of folders to run Openfoam;
  2. Export the created mesh to unv;
  3. Convert the unv to Openfoam;
  4. Prepare Icofoam to run.

1 - Prepare a set of folders to run Openfoam

I'm here not to reinvent the wheel. You can see a good Openfoam introduction in http://www.openfoam.org/docs/user/cavity.php

At this type of annular incompressible flow I choose to use Icofoam ... it is not the closest flow model to the problem ... it is the flow model. So, pick the cavity problem itself and copy to where you are going to run your simulation.

2 - Export the created mesh to unv

The best way, to me, is to export to a unv file. 
Right click over the mesh and choose to export to UNV file
When you export the file, please, export to the root folder of your problem. For example, if you copied "cavity" folder, this mesh should be on "cavity".

3 - Convert unv file to Openfoam

If you read some of my reports, you will see, we set the environment to succeed on this process. We have the Inlet, Outlet, Inner Wall and Outer Wall boundaries. So, this process should be straight forward. If you don't succeed you are maybe facing issues in your mesh or in your Salomé (6.5!) version (Take a look at http://openfoamwiki.net/index.php/IdeasUnvToFoam).

Is your unv file in the correct folder? If yes, just type:

  ideasUnvToFoam Annular.unv

My file is called "Annular.unv"

Ok. Type:
  
  checkMesh

And in your output you should be looking for this:

Checking patch topology for multiply connected surfaces ...
    Patch                      Faces    Points   Surface topology                  
    Inlet                       1296      1368     ok (non-closed singly connected)  
    Outlet                    1296      1368     ok (non-closed singly connected)  
    Outer_Wall           14400    14472    ok (non-closed singly connected)  
    Inner_Wall            14400    14472    ok (non-closed singly connected)  

As you can see, we have everything. Inlet, Outlet, Outer_Wall and Inner_Wall.

4 - Prepare Icofoam to run

Now comes the part that took me a while to understand. Take a look on the following link http://www.openfoam.org/docs/user/boundaries.php

Let's take care of the boundary conditions. Open the "p" file on the subfolder "0" (zero) and substitute for this one:


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

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

internalField   uniform 0;

boundaryField
{
    Inlet
    {
        type            zeroGradient;
    }

    Outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    Inner_Wall
    {
        type            zeroGradient;
    }

    Outer_Wall
    {
        type            zeroGradient;
    }

}

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

You can find an explanation of this file in here: http://www.openfoam.org/docs/user/basic-file-format.php.

But the there are some few things that I must call your attention.

On the line that says:

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

Is saying that you have a dimension of m²/s². Well, it is a pressure file so, we should be expecting in SI units a dimension of kg/(m.s²). What is the problem? No problem at all. With Icofoam we will use kinematic viscosity and to keep the units coherent on the solver we must divide the pressure by density.

Type zeroGradient:

As you can see, I'm using zeroGradient in every patch except "Outlet". It is only because I'm not using the gravity for this problem and solving the Navier-Stokes Equations for this setting you will find a zero pressure gradient. On the Outlet I'm using fixed value because the flow runs is not constrained anymore and goes to the environment that has zero pressure.

On the same "0" subfolder let us see the "U" file and substitute for the following:

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

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

internalField   uniform (0 0 0);

boundaryField
{
    Inlet
    {
        type            fixedValue;
        value           uniform (0 0 1.66);
    }

    Outlet
    {
        type            zeroGradient;
    }
 
    Inner_Wall
    {
type            fixedValue;
        value           uniform (0 0 0);
    }

    Outer_Wall
    {
type            fixedValue;
        value           uniform (0 0 0);
    }
}

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


As you can see, is very straightforward. You have the no-slip condition on the walls and the velocity inlet.

Now we should move to the subfolder "constant". In there you should find the polyMesh folder (don't need to touch it. It came from the unv conversion) and "transportProperties" file. Open it and substitute for:

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

nu              nu [ 0 2 -1 0 0 0 0 ] 7.27e-5;


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

As I said, IcoFoam uses the kinematic viscosity. You have the correct dimensions and the attribute value that came from http://chasingaftermystuff.blogspot.com.br/2013/05/using-icofoam-to-solve-annular-flow.html.

To the folder "System" I would advise you to wait a little bit more. I don't understand it as a whole yet and would suggest you to look into Openfoam's website for further information.

But if you have any doubt, just contact me ...




quinta-feira, 6 de junho de 2013

Creating a Structured Annular Grid with Viscous Layers in Salomé

Hi, from the previous posts you should be able to read the following python code to create your structured annular grid.

What is different from previous posts? Right now the code is very simple and we have viscous layers on the walls.

I had learned how to do this by taking a look at Salomé's forums. But, the function that called my attention and made this possible was PROPAGATE (http://docs.salome-platform.org/salome_7_2_0/gui/GEOM/geompy_doc/group__l3__blocks__op.html#ga81185dc6c66632dff79c2f0a19f77537)

With this function I could divide my inlet/outlet in a way that I could use the quadrangle mesh.

Again, if you have any doubt, just send me an email.


Structured Annular Mesh with Viscous Layer


###########
import sys
import salome

salome.salome_init()
theStudy = salome.myStudy

import salome_notebook
notebook = salome_notebook.notebook

###
### GEOM component
###

import GEOM
import geompy
import math
import SALOMEDS

gg = salome.ImportComponentGUI("GEOM")
geompy.init_geom(theStudy)

O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)

Wellbore = 8.5 # in
OD_Component = 6.5 # in
Length_Component = 10 # m

Wellbore = geompy.MakeCylinder(O, OZ, 0.5*Wellbore, Length_Component)
OD_Component = geompy.MakeCylinder(O, OZ, 0.5*OD_Component, Length_Component)
Annular_1 = geompy.MakeCut(Wellbore, OD_Component)
Annular_2 = geompy.MakeScaleAlongAxes(Annular_1, O, 0.0254, 0.0254, 1) # converting OD to imperial units

# Preparing for mesh
Plane_1 = geompy.MakePlaneLCS(None, 2000, 3)
Annular = geompy.MakeHalfPartition(Annular_2, Plane_1)
[Compound_1, Compound_2, Compound_3, Compound_4] = geompy.Propagate(Annular)
[Face_1,Face_2,Face_3,Face_4,Face_5,Face_6,Face_7,Face_8,Face_9,Face_10] = geompy.ExtractShapes(Annular, geompy.ShapeType["FACE"], True)
listSubShapeIDs = geompy.SubShapeAllIDs(Annular, geompy.ShapeType["FACE"])
Inlet = geompy.CreateGroup(Annular, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet, [42, 26])
Outlet = geompy.CreateGroup(Annular, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet, [45, 14])
Outer_Wall = geompy.CreateGroup(Annular, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outer_Wall, [4, 38])
Inner_Wall = geompy.CreateGroup(Annular, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inner_Wall, [48, 34])
Inlet_Outlet_Wire = geompy.CreateGroup(Annular, geompy.ShapeType["COMPOUND"])
geompy.UnionList(Inlet_Outlet_Wire, [Compound_3, Compound_4])

geompy.addToStudy( Wellbore, 'Wellbore' )
geompy.addToStudy( OD_Component, 'OD_Component' )
geompy.addToStudy( Annular, 'Annular' )
geompy.addToStudyInFather( Annular, Compound_1, 'Compound_1' )
geompy.addToStudyInFather( Annular, Compound_2, 'Compound_2' )
geompy.addToStudyInFather( Annular, Compound_3, 'Compound_3' )
geompy.addToStudyInFather( Annular, Compound_4, 'Compound_4' )
geompy.addToStudyInFather( Annular, Inlet, 'Inlet' )
geompy.addToStudyInFather( Annular, Outlet, 'Outlet' )
geompy.addToStudyInFather( Annular, Outer_Wall, 'Outer_Wall' )
geompy.addToStudyInFather( Annular, Inner_Wall, 'Inner_Wall' )


###
### SMESH component
###

import smesh, SMESH, SALOMEDS

smesh.SetCurrentStudy(theStudy)
import StdMeshers

Mesh_1 = smesh.Mesh(Annular)

Z_axis_segments = 200
Radius_segments = 10
Radius_divisions = 36 # For each half of inlet/outlet - It is going to divide 180 degrees for Radius_divisions

# Number of segments in z axis
Regular_1D = Mesh_1.Segment()
Nb_Segments_1 = Regular_1D.NumberOfSegments(Z_axis_segments)
Nb_Segments_1.SetDistrType( 0 )

# Number of radius segments
Regular_1D_1 = Mesh_1.Segment(geom=Inlet)
Nb_Segments_2 = Regular_1D_1.NumberOfSegments(Radius_segments)
Nb_Segments_2.SetDistrType( 0 )
Regular_1D_2 = Mesh_1.Segment(geom=Outlet)
status = Mesh_1.AddHypothesis(Nb_Segments_2,Outlet)

# Number of radius divisions
Regular_1D_3 = Mesh_1.Segment(geom=Inlet_Outlet_Wire)
Nb_Segments_3 = Regular_1D_3.NumberOfSegments(Radius_divisions)
Nb_Segments_3.SetDistrType( 0 )

# Type of meshing - Quadrangle for this case
Quadrangle_2D = Mesh_1.Quadrangle(algo=smesh.QUADRANGLE)
Quadrangle_Parameters_1 = Quadrangle_2D.QuadrangleParameters(StdMeshers.QUAD_STANDARD)

Quadrangle_2D_1 = Mesh_1.Quadrangle(algo=smesh.QUADRANGLE,geom=Inlet)
status = Mesh_1.AddHypothesis(Quadrangle_Parameters_1,Inlet)

Quadrangle_2D_2 = Mesh_1.Quadrangle(algo=smesh.QUADRANGLE,geom=Outlet)
status = Mesh_1.AddHypothesis(Quadrangle_Parameters_1,Outlet)

# Creation of a hexahedron meshing - it is going to extrude the base along z axis
Hexa_3D = Mesh_1.Hexahedron(algo=smesh.Hexa)
Viscous_Layers_1 = Hexa_3D.ViscousLayers(0.001,4,1,[ 42, 26, 45, 14 ]) # Viscous layer

isDone = Mesh_1.Compute()
Inlet_1 = Mesh_1.GroupOnGeom(Inlet,'Inlet',SMESH.FACE)
Outlet_1 = Mesh_1.GroupOnGeom(Outlet,'Outlet',SMESH.FACE)
Outer_Wall_1 = Mesh_1.GroupOnGeom(Outer_Wall,'Outer_Wall',SMESH.FACE)
Inner_Wall_1 = Mesh_1.GroupOnGeom(Inner_Wall,'Inner_Wall',SMESH.FACE)

## set object names
smesh.SetName(Mesh_1.GetMesh(), 'Mesh_1')
smesh.SetName(Regular_1D.GetAlgorithm(), 'Regular_1D')
smesh.SetName(Nb_Segments_1, 'Nb. Segments_1')
smesh.SetName(Quadrangle_2D.GetAlgorithm(), 'Quadrangle_2D')
smesh.SetName(Quadrangle_Parameters_1, 'Quadrangle Parameters_1')
smesh.SetName(Hexa_3D.GetAlgorithm(), 'Hexa_3D')
smesh.SetName(Viscous_Layers_1, 'Viscous Layers_1')
smesh.SetName(Nb_Segments_2, 'Nb. Segments_2')
smesh.SetName(Nb_Segments_3, 'Nb. Segments_3')
smesh.SetName(Inlet_1, 'Inlet')
smesh.SetName(Outlet_1, 'Outlet')
smesh.SetName(Outer_Wall_1, 'Outer_Wall')
smesh.SetName(Inner_Wall_1, 'Inner_Wall')

if salome.sg.hasDesktop():
  salome.sg.updateObjBrowser(1)
##############

quarta-feira, 5 de junho de 2013

Back to Business

Ufa ... this last embark was quite something. We had two influxes: one from water and another from gas. It has been some time that I didn't participate of a kick. But, thanks for the people we had on the rig (Toolpusher, CoMan, Drillers and us, DD's) we had identified everything from the beginning. A real team work. No harms were done to the people nor to the environment.

There are some few things I observed that I would like to test using CFD:
  • How the pump pressure is going to behave when we have a lighter mud on the annular?
  • How the pump pressure is going to behave when we have a heavier mud on the annular?
  • How the gas is going to behave inside of the annular while migrating?
  • How mud weight acts on gas migration? 
I didn't have time to search on SPE papers these items, and they do exist in there ... but ... I want to try to solve them using CFD.

While on the rig I couldn't run any simulation but I had studied Salomé a little bit more. Now I'm able to create a quite simple annular structured mesh with viscous layers. 

So, to a next post I'll use this method and solve again the annular problem using Openfoam. Expect to find on the next two posts:
  • Creating a Structured Annular Grid with Viscous Layers in Salomé;
  • How to set Icofoam to solve the Annular Flow problem.

Visit from the sky.
Dolphins trying to hijack the rig. We had to use water hoses to disperse them. :P