[DoQuocVu]

Hi everyone,

I'm trying to simulate a channel flow with cyclic boundary condition.
First, i was able to add an extra vector to the momentum equation in order to keep a constant pressure gradient along the channel.
Now i want to keep the cyclic boundary condition and add a cylinder inside the channel to see how the vortex distribute and i was trying to define the cylinder by using Immersed Boundary Method. So i made a few adjustment in the tutorial cylinderInChannelIcoIbFoam which is available in foam-extend 3.2.

However i have faced this error message when running icoIbFoam.

Create time

Create mesh for time = 0

Create immersed boundary cell mask
Create immersed boundary face mask
Found immersed boundary patch 0 named ibCylinder
External flow
Number of IB cells: 39
Reading field p

Reading field U


Calculating potential flow


--> FOAM FATAL ERROR: 
Can't find nearest triSurface point for cell 1800, (-0.98 0.48 0.05)Hit data = 0 (0.0528086 0.243847 0.05) -1


    From function immersedBoundaryFvPatch::makeIbPointsAndNormals() const
    in file immersedBoundaryFvPatch/immersedBoundaryFvPatch.C at line 1035.

FOAM aborting

Aborted (core dumped)

i was going through the immersedBoundaryFvPatch.C to see how the code working and my hypothesis is that the cyclic boundary condition made impact on the ibCellCentres[] list so that by some reason the code cannot find the cell that hit the triSurface. However it is still unclear because of my limited understanding in IBM method and programming.

i would very appreciate any of your help to solve this error.

[hjasak]

You probably have a piece of STL boundary outside of your mesh.

Hrv

[DoQuocVu]

Dear prof Jasak,

Thank you for your reply.
However, all i did is to change the boundary condition of inlet and outlet to cyclic type and did not adjust the STL mesh.

I've read through the immersedBoundaryFvPatch.C and printed out all the ibCells to check. By that way, I figured out within the list there were two invalid cells at the top-left corner and top-right corner, since they are quite far from the cylinder surface (named 1800 and 1874). Therefore, i guess the problem is probably the cyclic boundary condition that made the code confuse when checking immersed boundary cells.

I also tried hiding these following lines of the immersedBoundaryFvPatch.C (from line 420 to 467) and printed out ibCells again and this time the two cells above were gone.

    forAll (gE.boundaryField(), patchI)
    {
        if (gE.boundaryField()[patchI].coupled())
        {
            scalarField gammaExtOwn =
                gE.boundaryField()[patchI].patchInternalField();

            scalarField gammaExtNei =
                gE.boundaryField()[patchI].patchNeighbourField();

            const unallocLabelList& fCells =
                mesh_.boundary()[patchI].faceCells();

            forAll (gammaExtOwn, faceI)
            {
                if
                (
                    mag(gammaExtNei[faceI] - gammaExtOwn[faceI])
                  > SMALL
                )
                {
                    if (gammaExtOwn[faceI] > SMALL)
                    {
                        if (!ibCellSet.found(fCells[faceI]))
                        {
                            ibCellSet.insert(fCells[faceI]);
                        }
                    }
                    else if (2*gammaExtOwn.size() == fCells.size())
                    {
                        if
                        (
                           !ibCellSet.found
                            (
                                fCells[gammaExtOwn.size() + faceI]
                            )
                        )
                        {
                            ibCellSet.insert
                            (
                                fCells[gammaExtOwn.size() + faceI]
                            );
                        }
                    }
                }
            }
        }
    }

Furthermore the origin tutorial still can run smoothly without the above lines of code. However, i'm not fully understand how the code work and so i have no idea what to do now.
So professor, could you give me some advice?

Vu

[DoQuocVu]

Dear prof Jasak,

It's been a while and i'm still struggling with applying cyclic boundary condition to immersed boundary solver. I have some ideas of how the code works now.

In the makeIbCells() function, there are 2 forAll loops. The first one is used to check interior ibCells and the second is to check on boundary. If the boundary is coupled, the code will take gamma values of patchInternalField and patchNeighbourField into account.

The problem is the gamma values returned by gE.boundaryField().patchNeighbourField() is always wrong. Specifically, In the cylinderInChannelIcoIbFoam tutorials, the cylinder is totally located within internalField, so the gamma values at centroids of all boundaries's cell must be equal to 1. However, at the second forAll loops, i printed out the gamma values which is assigned to gammaExtOwn and gammaExtNei:

forAll (gE.boundaryField(), patchI)
    {
        if (gE.boundaryField()[patchI].coupled())
        {    
              scalarField gammaExtOwn =gE.boundaryField()[patchI].patchInternalField();
              Pout<< "-----gamma values at owner cells: " << gammaExtOwn << endl; // Added
              scalarField gammaExtNei =gE.boundaryField()[patchI].patchNeighbourField();
              Pout<< "-----gamma values at neighbour cells: " << gammaExtNei << endl; // Added

              const unallocLabelList& fCells =
                        mesh_.boundary()[patchI].faceCells();
...

Here is what i got:

-----gamma values at owner cells: 25{1}
-----gamma values at neighbour cells: 
25
(
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0.004
)

-----gamma values at owner cells: 25{1}
-----gamma values at neighbour cells: 
25
(
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0
)

The 0.004 and 0 values returned by patchNeighbourField is clearly unreasonable in this situation. I would like to find the implementations of patchNeighbourField to figure out what exactly it does but i couldn't. So i stuck again. Do you have any idea what is happening here inside the code?

[Ben UWIHANGANYE]

Hello Foamers,

May i get your help? Indeed, I need to simulate 2-D oscillating cylinder with immersed boundary Mesh. Up to date, I have been able to write simple harmonic equation of motion in C++ (x''+w^2x=0). But the challenging task now is to map a 2-d cylinder within Open Foam that will be oscillating and
again to mesh it is with IB not body fitted meshes.

Can you help to know through which way, I will be able to make that solver?

Regard!

[DoQuocVu]

Hi Ben,

Have you gone through the movingCylinderInChannelIco tutorial? I think to make the immersed cylinder oscillate, you only need to define the motion in dynamicMeshDict the same way as you do with a moving body fitted cylinder

[Ben UWIHANGANYE]

Originally Posted by DoQuocVu

Hi Ben,

Have you gone through the movingCylinderInChannelIco tutorial? I think to make the immersed cylinder oscillate, you only need to define the motion in dynamicMeshDict the same way as you do with a moving body fitted cylinder

It doesn't work for me. I have tried to install foam extend 4.0 but failed.

There is no other way I can edit an existing solver like pimpleDyMFoam and get immersed boundary meshes. Just guide me!

Regard!

[Ben UWIHANGANYE]

Hello,

I wanted to plot forces coefficients for a movingcylinderInChannel tutorial. Can you help?

how should my controldict look like?

I would appreciate your help.

Thank you!

[DoQuocVu]

You can simply use functionObjects, add these lines to the end of controlDict
functions
(
forces
{
type immersedBoundaryForces;
functionObjectLibs ("libimmersedBoundaryForceFunctionObject.so");
outputControl timeStep;
outputInterval 1;
patches ( nameOfyourpatch);
pName p;
UName U;
rhoName rhoInf;
rhoInf 1;
log true;
CofR ( 0 0 0 );
Aref 0.05; //Your cylinder's diameter if 2D and cylinder's surface area if 3D
Uref 1; // Inlet velocity
}
);
just take a look at the controlDict of the movingCylinderInChannel tutorials and do the exact same thing!

Hope this can help

[Ben UWIHANGANYE]

Originally Posted by DoQuocVu

You can simply use functionObjects, add these lines to the end of controlDict
functions
(
forces
{
type immersedBoundaryForces;
functionObjectLibs ("libimmersedBoundaryForceFunctionObject.so");
outputControl timeStep;
outputInterval 1;
patches ( nameOfyourpatch);
pName p;
UName U;
rhoName rhoInf;
rhoInf 1;
log true;
CofR ( 0 0 0 );
Aref 0.05; //Your cylinder's diameter if 2D and cylinder's surface area if 3D
Uref 1; // Inlet velocity
}
);
just take a look at the controlDict of the movingCylinderInChannel tutorials and do the exact same thing!

Hope this can help

Dear Vu,

I thank you for your quick feedback.
However, i wanted the coefficients(drag coefficients and lift coeff.). Can you help on that? or it is done in the same way?

Another point, i want to simulate similar cylinder (D=0.05m) but at different Reynolds number. Can you help to know the Reynold used in that tutorial and how i can modify it.

I would appreciate.

Regard!

[Ben UWIHANGANYE]

Hello foamers,

I was running the tutorial in Foam extent 4.0(movingcylinderInchannel) but failed to refine mesh in the vicinity of the cylinder.

Can any one help to know how I can refine the mesh.


Thank you!

[Ben UWIHANGANYE]

Hello

I would like to ask about convergence of a simulation in foam ext 4.0
Indeed, i am running a case following a tutorial of immersed boundary(movingcylinderInchannel), using IcoDyMIbFOAM,

I would like to ask about its convergence, Does it converge, according to you what should I do to get it converged?

You may refer to the attached picture of residuals of my case!


Regard!
.................................................. ......................................
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
p
{
solver amgSolver;
cycle W-cycle;
policy PAMG;
nPreSweeps 2;
nPostSweeps 2;
groupSize 4;
minCoarseEqns 4;
nMaxLevels 100;
scale on;
smoother symGaussSeidel;

minIter 1;
maxIter 100;
tolerance 1e-7;
relTol 0.01;
}

pFinal
{
solver amgSolver;
cycle W-cycle;
policy PAMG;
nPreSweeps 2;
nPostSweeps 2;
groupSize 4;
minCoarseEqns 4;
nMaxLevels 100;
scale on;
smoother symGaussSeidel;

minIter 1;
maxIter 100;
tolerance 1e-7;
relTol 0.00;
}

U
{
solver BiCGStab;
preconditioner ILU0;

minIter 1;
maxIter 1000;
tolerance 1e-08;
relTol 0;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 5;

pRefPoint (0 -0.45 0.05);
pRefValue 0;
}

PIMPLE
{
nOuterCorrectors 1;
nCorrectors 4;
nNonOrthogonalCorrectors 0;

pRefPoint (0 -0.45 0.05);
pRefValue 0;
}

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


and

FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application icoDyMIbFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 40;
deltaT 0.02;
writeControl adjustableRunTime;
writeInterval 0.5;
purgeWrite 0;
writeFormat binary;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
adjustTimeStep yes;
maxCo 0.5;
libs
(
"libimmersedBoundary.so"
"libimmersedBoundaryDynamicFvMesh.so"
);
functions
(
forces
{
type immersedBoundaryForces;
functionObjectLibs ("libimmersedBoundaryForceFunctionObject.so");

outputControl timeStep;
outputInterval 1;
patches ( ibCylinder );

pName p;
UName U;
rhoName rhoInf;
rhoInf 1;

log true;
CofR ( 0 0 0 );

Aref 0.05;
Uref 20;
}

);


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