CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Improve simpleFoam convergence

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree6Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   July 22, 2010, 09:52
Default
  #21
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Quote:
Dear Matey
But if i inizialized my simulation with a field very precise, I must not use necessarily a very low tolerance or not?
Yes, you should do! The precise field for initialization gets "unprecise" immediately on the new mesh.

Quote:
Post my solution. The source for mapfields utility is the same for all. I change tolerance convergence.
Although your results for tau/tau_0 look stable after reaching 1e-2, it really is a very high tolerance criteria! I use at least 1e-4 for simpleFoam. The other results (U for example) will differ much more, I suppose. Well, I hope you are not running nuclear power plants or passenger air planes with your simulations
Alhasan likes this.
MartinB is offline   Reply With Quote

Old   July 22, 2010, 10:58
Default
  #22
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
No I'm not simulating a nuclear power plant! But also u profile is good, i think that I'll use 1e-3 tolerance convergence. I know the importance of tolerance i made structural simulation for a wing with nastran and I know the its importance
Daniele111 is offline   Reply With Quote

Old   July 22, 2010, 12:15
Default
  #23
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
Hi
Do you know a way to interface my openfoam simulation with a matlab script that use OpenFoam results?
Daniele111 is offline   Reply With Quote

Old   July 22, 2010, 14:46
Default
  #24
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
The best interface is interaction via files, I think...
I found this thread where you had problems starting OpenFOAM from Matlab:
OpenFOAM command from inside MATLAB

Do you still have problems in doing so? If so, try
Code:
    unix('foamExec simpleFoam')
since this should set your environment variables properly.

Otherwise
Code:
    unix('source path/to/OpenFOAMs/bashrc; simpleFoam')
could work.

For the data you need in Matlab:
Do you only need values for shear stress on the wall?
Is it enough to have coordinates for each shear stress value, or do you need the full mesh structure?
Do you use a more sophisticated transport model than newtonian flow?
Are you using OpenFOAM in parallel?

I have some code snippets here that would write a simple ";"-separated .csv file which you can parse in Matlab.
MartinB is offline   Reply With Quote

Old   July 22, 2010, 15:08
Default
  #25
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
If I use your command I have this:

unix('foamExec blockMesh')
/bin/bash: foamExec: command not found


I only need values for shear stress on the wall, and it is enough to have coordinates for each shear stress value;

"Do you use a more sophisticated transport model than newtonian flow?" ->> No
Are you using OpenFOAM in parallel? ->> No (maybe later)

Thanks
Daniele
Daniele111 is offline   Reply With Quote

Old   July 22, 2010, 15:24
Default
  #26
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
What about the other command: unix('source path/to/OpenFOAMs/bashrc; simpleFoam') ?

"path/to/OpenFOAMs" must be replaced with the full path to your OpenFOAM's installation directory, than the "bashrc" file follows, which includes the environment for OpenFOAM. Than there is a ";", followed by another command (simpleFoam) to be executed in the just set environment...

What's the name of the patch, you need the shear stress values for? Is it "bottom"?
MartinB is offline   Reply With Quote

Old   July 22, 2010, 15:33
Default
  #27
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
I do:
unix('/home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc; blockMesh')
and

/bin/bash: /home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc: permission denied
/bin/bash: blockMesh: command not found

Patch name is "bottom"

How long do you work with OpenFoam? Do you use it for university or work? I'm star to use is one month ago, and I like it, it's very flexible.
Daniele111 is offline   Reply With Quote

Old   July 22, 2010, 16:37
Default
  #28
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
unix('/home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc; blockMesh')

Now It works, I changed matlab permission

Thanks
Daniele111 is offline   Reply With Quote

Old   July 22, 2010, 16:50
Default
  #29
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Daniele,

here we go:

Go to the simpleFoam source directory:
"~/OpenFOAM/OpenFOAM-1.6.x/applications/solvers/incompressible/simpleFoam/"
The version number might be different for your system.

Create a new file there with name "write.H". Paste this code in there:
Code:
// write matlab file
{
    Info << "Writing wall shear stress for matlab!" << endl;

    word wallPatchName = "bottom";
    Info << "Searching for patch " << wallPatchName << endl;

    // Check, if the patch name exists:
    label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
    if (patchID < 0)
    {
        Info << "No patch found with name \"" << wallPatchName << "\"!"
             << endl;
    }
    else
    {
        /* for laminar case:
        Info<< "Reading viscosity nu\n" << endl;
        dimensionedScalar nu_
        (
            laminarTransport.lookup("nu")
        );
        */
        // calculating wall shear rate
        volVectorField wallGradU
        (
                IOobject
                (
                        "wallGradU",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedVector
                (
                        "wallGradU",
                        U.dimensions()/dimLength,
                        vector::zero
                )
        );

        forAll(wallGradU.boundaryField(), patchi)
        {
            wallGradU.boundaryField()[patchi] =
                -U.boundaryField()[patchi].snGrad();
        }

        // calculating wall shear stress
        volVectorField tau
        (
                IOobject
                (
                        "tau",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                //laminar case: nu_ * wallGradU
                turbulence->nuEff() * wallGradU
                // turbulence->nu() * wallGradU
        );

        // Pointer to the patch we want to have shear stress for:
        const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
        label nFaces = cPatch.size();

        pointField faceCenterCoords(nFaces, vector(0,0,0));
        pointField faceAreas(nFaces, vector(0,0,0));
        pointField wallShearStress(nFaces, vector(0,0,0));

        label runFaces = 0;

        // collecting the interesting data
        for (int i=0; i<nFaces; i++)
        {
            faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
            faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
            wallShearStress[runFaces] = wallGradU.boundaryField()[patchID][i];
            runFaces++;
        }

        //  make a cute filename
        fileName casePath = cwd();//casePath.name();
        fileName caseName = casePath.name();
        fileName myFile = caseName+"_wallShearStress.csv";
        Info << "Data is written to file " << myFile << endl;

        // create a new file
        OFstream *myStream = NULL;
        myStream = new OFstream(myFile);

        for (int i=0; i<nFaces; i++)
        {
            *myStream << i << "; "
                    // coordinates on face centers:
                    << mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
                    << mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
                    << mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

                    // wall shear stress as a vector:
                    << tau.boundaryField()[patchID][i].component(0) << "; "
                    << tau.boundaryField()[patchID][i].component(1) << "; "
                    << tau.boundaryField()[patchID][i].component(2) << "; "

                    // magnitude of wall shear stress:
                    << mag(tau.boundaryField()[patchID][i]) << endl;
        }

    }
}
Then open the simpleFoam.C file. Add this line to the includes at the top:
Code:
#include "OFstream.H"
Behind the runTime loop of simpleFoam paste the following two lines, so that it looks like:
Code:
        ...
        #include "convergenceCheck.H"
    }
    // writing additional data for matlab
#   include "write.H" // <--- this is the important one
    ...
In the simpleFoam directory call
Code:
wclean
wmake
Now if you are running simpleFoam an extra file with name "your_case_name_wallShearStress.csv" is written. It's format is:
number; coordX; coordY; coordZ; wallShearStressX; wallShearStressY; wallShearStressZ; magnitude(wallShearStress)

I haven't tested this code, so use it with care! Check and double check, if the wall shear stress is calculated really correctly!!!
Have a look a line 62, where the wall shear stress is calculated:
Code:
turbulence->nuEff() * wallGradU
is just a guess, since I don't have experience with turbulence. It might be nu, or nut or something totally different there...

Have fun :-)

Martin

p.s.: I'm using OpenFOAM for seven month now... but I'm still learning, of course ;-) And I'm using it for work.
MartinB is offline   Reply With Quote

Old   July 23, 2010, 04:04
Default
  #30
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
Dear Martin
Thanks for your code, is a gold mining! There are many things I want to learn to. I'll change the computation of wallShearStress, I think that i'll use the openfoam code for this and I add it to write.h.

Daniele
Daniele111 is offline   Reply With Quote

Old   July 26, 2010, 11:05
Default MRF simple foam linear upwind: convergence troubles
  #31
Member
 
Aldo Iannetti
Join Date: Feb 2010
Posts: 48
Rep Power: 7
aldo.iannetti is on a distinguished road
Hi,
I'm trying to set a linearUpwind MRF simpleFoam simulation of a cooling fan. The mesh has been imported from GAMBIT (8.5 e+6 cells, not bad quality) and the simulation diverges.Here attached: log, fvSchemes,fvSolution. Could you please have a look in order to give me advices about settings?
thanks
Aldo
Attached Files
File Type: zip MRFSImpleFoam_linearUpwind.zip (96.1 KB, 19 views)
aldo.iannetti is offline   Reply With Quote

Old   August 21, 2010, 06:24
Default
  #32
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
But if I want extract also k with code? How can i modify the code post previously?
Daniele111 is offline   Reply With Quote

Old   August 21, 2010, 08:37
Default
  #33
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Daniele,

what about adding
Code:
<< k.boundaryField()[patchID][i] << "; "
to the "myStream" loop?
Martin
MartinB is offline   Reply With Quote

Old   August 21, 2010, 08:47
Default
  #34
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
Hi Martin
It is that I do...I did a syntax error probably.
Thanks!

Daniele
Daniele111 is offline   Reply With Quote

Old   November 4, 2010, 20:12
Default
  #35
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
Hi
Now I try to use a cluster to increase my speed, and so i'm traying tu use OpenFoam parallel utility. But when i run decomposePar i have this error:

FOAM FATAL IO ERROR:

Cannot find 'value' entry on patch bottom of field epsilon in file "/u/acconcia/OpenFOAM/acconcia-1.6/run/blayer4/0/epsilon"
which is required to set the values of the generic patch field.
(Actual type epsilonWallFunction)

Please add the 'value' entry to the write function of the user-defined boundary-condition
or link the boundary-condition into libfoamUtil.so

file: /u/acconcia/OpenFOAM/acconcia-1.6/run/blayer4/0/epsilon::boundaryField::bottom from line 49 to line 52.

From function genericFvPatchField<Type>::genericFvPatchField(con st fvPatch&, const Field<Type>&, const dictionary&)
in file genericFvPatchField/genericFvPatchField.C at line 71.


My boundary condition are the same of the case that I use not in parallel and they haven't give me problem.
Daniele111 is offline   Reply With Quote

Old   November 5, 2010, 00:50
Default
  #36
Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 92
Rep Power: 8
matejfor is on a distinguished road
The best think you can do is to follow what the error message clearly says so (especially when it does say it so clearly). On patch "bottom" you have not prescribed value of epsilon. What is the BC for epsilon on bottom? It is probably some which needs a value.

Now why you do not need it in single but need it in parallel?
My guess is that you do not use k-eps but k-omega where BC for epsilon is not needed.
The decompose utility does not know what you will use for the run so it does decompose everything and checks for you BCs which is very kind of it.
But it's only my guess.


good luck

matej
matejfor is offline   Reply With Quote

Old   November 5, 2010, 04:45
Default
  #37
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
These are my BC. And I'm sure that I'm using k-eps model, this BC have never been a problem.

Thanks
Attached Files
File Type: gz 0.tar.gz (1.3 KB, 19 views)
Daniele111 is offline   Reply With Quote

Old   November 5, 2010, 05:04
Default
  #38
Member
 
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 92
Rep Power: 8
matejfor is on a distinguished road
OK. then you have 2 possibilities:

(1) add the value uniform number; line to your bottom BC for the decomposer to be happy
(end remove it afterwards as the value is not needed by the wall functions (you can check it in doxygen)

(2) remove the need of it from the decomposer, recompile it and if it works report it


matej
matejfor is offline   Reply With Quote

Old   November 5, 2010, 11:40
Default
  #39
Senior Member
 
Daniele
Join Date: Feb 2010
Posts: 134
Rep Power: 7
Daniele111 is on a distinguished road
Hi
I solved my problem, paralle is happy with a value!
To calculate wall shear stress with single processor I used Martin's code of post 29, but in parallel I can't use it: is possible a similar code to have a output at simulation stop? A tau.dat like output then i'll manipulate with matlab.
Daniele111 is offline   Reply With Quote

Old   November 5, 2010, 15:16
Default
  #40
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 252
Rep Power: 11
MartinB is on a distinguished road
Hi Daniele,

not fully tested, but it should be a point to start from:

Code:
// write matlab file
{
    Info << "Writing wall shear stress for matlab!" << endl;

    word wallPatchName = "bottom";
    Info << "Searching for patch " << wallPatchName << endl;

    // Check, if the patch name exists:
    label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
    if (patchID < 0)
    {
        Info << "No patch found with name \"" << wallPatchName << "\"!"
             << endl;
    }
    else
    {
        /* for laminar case:
        Info<< "Reading viscosity nu\n" << endl;
        dimensionedScalar nu_
        (
            laminarTransport.lookup("nu")
        );
        */
        // calculating wall shear rate
        volVectorField wallGradU
        (
                IOobject
                (
                        "wallGradU",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedVector
                (
                        "wallGradU",
                        U.dimensions()/dimLength,
                        vector::zero
                )
        );

        forAll(wallGradU.boundaryField(), patchi)
        {
            wallGradU.boundaryField()[patchi] =
                -U.boundaryField()[patchi].snGrad();
        }

        // calculating wall shear stress
        volVectorField tau
        (
                IOobject
                (
                        "tau",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                //laminar case: nu_ * wallGradU
                turbulence->nuEff() * wallGradU
                // turbulence->nu() * wallGradU
        );
        
        // Pointer to the patch we want to have shear stress for:
        const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
        label nFaces = cPatch.size();

        pointField faceCenterCoords(nFaces, vector(0,0,0));
        pointField faceAreas(nFaces, vector(0,0,0));
        pointField wallShearStress(nFaces, vector(0,0,0));

        label runFaces = 0;

        // collecting the interesting data
        for (int i=0; i<nFaces; i++)
        {
            faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
            faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
            wallShearStress[runFaces] = tau.boundaryField()[patchID][i];
            runFaces++;
        }

        // for parallel
        reduce(nFaces, sumOp<label>());
        Info << "Number of global faces: " << nFaces << endl;
        

        //  make a cute filename
        fileName casePath = cwd();//casePath.name();
        fileName caseName = casePath.name();
        fileName myFile = caseName+"_wallShearStress.csv";
        Info << "Data is written to file " << myFile << endl;

        // create a new file
        OFstream *myStream = NULL;
        if (Pstream::master())
        {
            myStream = new OFstream(myFile);
        }

        // slaves send data to master, master writes data to file
        if (Pstream::parRun())
        {
            if (Pstream::myProcNo() != Pstream::masterNo())
            {
                OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                toMaster << faceCenterCoords;
                toMaster << faceAreas;
                toMaster << wallShearStress;
            }
            else // master part
            {
                // first write own data
                for (int i=0; i<faceCenterCoords.size(); i++)
                {
                    *myStream << i << ";"
                        << faceCenterCoords[i].component(0) << "; "
                        << faceCenterCoords[i].component(1) << "; "
                        << faceCenterCoords[i].component(2) << "; "

                        << wallShearStress[i].component(0) << "; "
                        << wallShearStress[i].component(1) << "; "
                        << wallShearStress[i].component(2) << "; "

                        << mag(wallShearStress[i]) << endl;
                }

                label run = faceCenterCoords.size();
                
                // then receive slave data and write
                pointField bufferFaceCenterCoords;
                pointField bufferFaceAreas;
                pointField bufferWallShearStress;

                
                for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                {
                    IPstream fromSlave(Pstream::blocking, slave);
                    fromSlave >> bufferFaceCenterCoords;
                    fromSlave >> bufferFaceAreas;
                    fromSlave >> bufferWallShearStress;
                    
                    for (int i=0; i<bufferFaceCenterCoords.size(); i++)
                    {
                        *myStream << run << ";"
                            << bufferFaceCenterCoords[i].component(0) << "; "
                            << bufferFaceCenterCoords[i].component(1) << "; "
                            << bufferFaceCenterCoords[i].component(2) << "; "

                            << bufferWallShearStress[i].component(0) << "; "
                            << bufferWallShearStress[i].component(1) << "; "
                            << bufferWallShearStress[i].component(2) << "; "

                            << mag(bufferWallShearStress[i]) << endl;
                        run++;
                    }
                }
            }
        }
        else // single processor only
        {
            for (int i=0; i<nFaces; i++)
            {
                *myStream << i << "; "
                        // coordinates on face centers:
                        << mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
                        << mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
                        << mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

                        // wall shear stress as a vector:
                        << tau.boundaryField()[patchID][i].component(0) << "; "
                        << tau.boundaryField()[patchID][i].component(1) << "; "
                        << tau.boundaryField()[patchID][i].component(2) << "; "

                        // magnitude of wall shear stress:
                        << mag(tau.boundaryField()[patchID][i]) << endl;
            }
        }

    }
}
Martin
MartinB is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Convergence Problems SimpleFOAM Kutti OpenFOAM 16 June 14, 2010 08:12
Getting faster convergence in simpleFoam basneb OpenFOAM 8 February 9, 2010 05:20
Definition of convergence criterion in simpleFoam titio OpenFOAM Running, Solving & CFD 1 February 6, 2010 02:34
Convergence of CFX field in FSI analysis nasdak CFX 2 June 29, 2009 01:17
to improve convergence Mohamed FLUENT 0 May 14, 2003 05:43


All times are GMT -4. The time now is 13:24.