[wenterodt]

Dear fellow users,
I am using simpleFoam for laminar calculations in axisymmetric domains. I set up the cases by the book, i.e.
-wedge shaped domain that straddles x-y-plane with an opening angle of 5deg
-only one cell in z-direction
-wedge BCs for front and back (separate patches)
-empty BC for axis.
The calculations converge beautifully, so far no problem.

BUT even though the axisymmetric problem is a 2d problem and Uz is initialized with 0, the solver tries to solve the system for the z-direction in every iteration. This does not only waste cpu-time, but it imposes an unnecessary error to the solution.

I would be very glad to know the answers to the following questions:
-Is this normal behavior for axisymmetric calculations or did I miss something in the setup?
-Can this be prevented in some way like using "empty" patches in plane 2d calculations?
-Are there naming conventions that have to be obeyed?
-Is it possible to manipulate the solver settings or the "solve(UEqn() == -fvc::grad(p))"-call to achieve this?

Thanks,
Tammo

[msrinath80]

I assume that you have already seen a related thread on this issue: http://www.cfd-online.com/OpenFOAM_D...tml?1211878793

[msrinath80]

Correct me if I am wrong, but isn't the Y-axis your third direction according to this mesh definition?

The following statement seems to indicate the same:

hex (0 1 2 3 0 1 4 5) (10 1 10) simpleGrading (1 1 1)

I don't have the time to go over your setup at this point, but I do get the feeling that you're doing something fundamentally wrong here.

[wenterodt]

Of course the blocks-line has to read like this:

blocks
(
hex (0 1 2 3 0 1 4 5) (10 10 1) simpleGrading (1 1 1)
);

But this does not solve the problem...

[shivasub]

Hi Tammo,

You can use a simple hack to make the Uz=0. This is how we use it in our code if we are running wedge cases.

// Remove the swirl component of velocity for "wedge" cases
if (piso.found("removeSwirl"))
{
label swirlCmpt(readLabel(piso.lookup("removeSwirl")));

U.field().replace(swirlCmpt, 0.0);
}

regards
Shiva

[wenterodt]

Dear Shiva,
thank you for your reply! I tried something alike before, but I wasn't aware that it is possible to set Uz=0 with only one line of code...

Setting Uz=0 in every iteration like you proposed does indeed help speeding up the iteration process since the number of iterations for Uz will decrease. However, I would like the solver not to bother about Uz at all since I am dealing with an axisymmetric (i.e. 2d) problem!

Does this hack prevent your solver from solving for Uz? In that case it would be clear that I set up my case wrongly, which would be good news, since that would propably be the easiest to correct...

As Hrvoje Jasak states here:
http://www.cfd-online.com/cgi-bin/Op...0697#POST10697
the "wedge"-BC will cause some transformations, but so far I can see no difference between using "wedge" and "cyclic" BC in an axisymmetric flow from the behavior of the solver.

I am puzzled...

Regards,
Tammo

[shivasub]

Hi Tammo,

yes, the hack provided from our code does not eliminate solving for Uz, it just reduces the number of iterations to something very small.

As Hrv says in his post, wedge is used only when the mesh is only one cell thick in the 3rd dimension, i.e the one face will map back to the opposite face of the same cell. Whilst cyclic will be used for periodic boundaries.

When you say you see no difference between using "wedge" and "cyclic", have you used both for axisymmetric cases and find no difference?

regards
Shiva

[wenterodt]

Shiva, I think my statement there was ambigous. What I meant to say was that I can not see why there must be "wedge" type boundary condition at all, at least how it is implemented in OpenFOAM.

From a physical point of view the cyclic bc should be sufficient to model swirl in axisymmetric flows, since all variables must be the same (i.e. can be mapped directly) on the front and on the back of the domain, especially Uz.

Before my first calculations with the wedge bc, I thought it would simplify the problem as to not solve for Uz, thus eliminate any swirl. This would seem reasonable to me (with 1-cell-thickness in z-direction):
-with swirl: cyclic
-without swirl: symmetryPlane
-without swirl (not solve for Uz): wedge

But (and this at last is the answer to your question) OpenFOAM behaves quite differently:
-wedge: models swirl
-symmetryPlane: no swirl, looks ok
-cyclic: gives results that are physically impossible

Still I am wondering if there is a way to model axisymmetric flow without swirl and eliminate the uneccessary solving for Uz in that case...

Regards,
Tammo

[benham]

Hi there,

I am struggling with the same issue: I can't get the "tangential" velocity component to converge when using the wedge boundary type.

It seems from above that the only way to fix it is by forcing the component to vanish at each iteration. (I don't know how to implement this.)

Is there really still no official fix for this?

Best regards,
Graham

[catelyn]

Hi everybody,

did anyone solved this issue?

Regards.

[benham]

No I never solved the issue. Whenever I use the wedge feature, I just run the solver for as many iterations it takes to get the other components to converge (Ux,Uy,k,epsilon,etc...). This is sort of annoying, but at least it is solving the problem correctly.

[francescomarra]

Dear All,

it looks to me that a basic misunderstanding avoids to properly consider the issue. A 2D field means only that, assuming it is 2D in the x,y plane, \partial \phi/\partial z = 0 for any variable \phi. Therefore Uz not equal zero and even varying with x and y are perfectly admissible in a 2D field and need to be resolved. This actually occurs in swirled flows.

I think that the code is correctly programmed to deal with this and that to avoid solving for Uz=0 has to be enforced in some special way, as for instance suggested by Shiva, when the user know from physics that its field has Uz=0 everywhere.

Best regards,

Franco

[bigphil]

Although not the most elegant, you can disable solving in a particular direction by changing the "solutionD()" vector, which is accessible (as a const reference) from the mesh: for example, you could add the following code to the start of your solver to disable solution in the z direction:

        Vector<label>& solD = const_cast<Vector<label>&>(mesh.solutionD());

        Info<< "Disabling solution in the z direction" << endl;
        solD[2] = -1;

Of course, you can add some code to determine the wedge direction automagically.

Philip

[benham]

Hi Philip,

Thanks for your help. Which file exactly should I put this code at the beginning?

I am using the simpleFoam solver on OF5.0.

Cheers,
Graham

[bigphil]

Originally Posted by benham

Hi Philip,

Thanks for your help. Which file exactly should I put this code at the beginning?

I am using the simpleFoam solver on OF5.0.

Cheers,
Graham

Hi Graham,

You can add this code just after '# include "createMesh.H"' in the simpleFoam.C file.

You then need to re-compile the solver (I suggest you make a copy of the solver and call it something like mySimpleFoam).

Philip

[Artur]

Hi All,


Thanks for the pointer Phillip, the modified version of the solver works like a charm, but I also agree with the earlier point by Franco that such arbitrary disabling of the solution direction is better left to the user since in certain cases it's not necessarily correct.


For those who are interested, here's a full example of how I modified the simpleFoam solver to control which direction should be disabled on an axisymmetric grid:

    #include "createMesh.H"
     // - Added code here. -
    int iNormalDir(-1);
    word normalDirection;
    
    IOdictionary transportPropertiesDict
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );
    transportPropertiesDict.lookup("normalDirection") >> normalDirection;
    
    if (normalDirection == "x")
        iNormalDir = 0;
    else if (normalDirection == "y")
        iNormalDir = 1;
    else if (normalDirection == "z")
        iNormalDir = 2;
    else
        FatalErrorIn(args.executable().c_str())
            << "Unknown Cartesian direction " << normalDirection << " attempting to be disabled!"
            << exit(FatalError);

    Vector<label>& solD = const_cast<Vector<label>&>(mesh.solutionD());
    Info << "Disabling solution in the " << normalDirection << " direction" << endl;
    solD[iNormalDir] = -1;
     // ---

I'm still running my benchmark case but it seems to be converging considerably faster judging by the residuals.


Happy foaming,


Artur

[joaran]

Hi all,


The code provided by Artur makes convergence a charm in just a few iterations. I noticed, though, that the solution is quite different than the one obtained by solving in all directions (and not achieving convergence in the normal velocity). When running a very simple case in OpenFOAM v9 using rhoPimpleFoam, where I add a small energy input through fvModels, I see that energy is conserved perfectly in the non-converged case, but when incorporating the suggested modification (which as I said before converges perfectly) the amount of energy in the domain is greatly increased (i.e. is not conserved at all).



Does anyone have a clue on why this happens?


Best,
Joaquín