CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Problems with simple test case (http://www.cfd-online.com/Forums/openfoam-solving/98835-problems-simple-test-case.html)

matthiass March 20, 2012 10:01

Problems with simple test case
 
3 Attachment(s)
Hello everybody,

I want to use OpenFoam for solving structural mechanic problems and at the moment I'm solving some test cases that I want to compare with analytical solutions.

My present case is very simple. It is a beam, which is clamped at one side. Then I put a surface load on the beam, so that it bends.
The analytical solution and the solution of a commercial FE-software are similar: displacement at the tip of about 1cm.

I attach my test case, so that someone who is more experienced than me with OpenFoam (which is not really difficult ;-)) can give me a hint.
The mesh of the *.unv is a bit coarse. But I already tried a mesh with 50k elements and got no better result

Best regards and thanks in advance for any reply
Matthias

P.S.: Here is my D-file, maybe there is already an error:
Code:

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

dimensions      [0 1 0 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{
    clamp
    {
        type            fixedValue;
    value          uniform (0 0 0);
    }
    move
    {
        type            tractionDisplacement;
        traction        uniform ( 0 10000 0) ;
        pressure        uniform 0;
    value            uniform (0 0 0);
    }
    free
    {
        type            tractionDisplacement;
        traction        uniform ( 0 0 0) ;
        pressure        uniform 0;
    value            uniform (0 0 0);
    }
}

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


kmooney March 20, 2012 10:48

Could you give a few more details on the loading configuration? i.e. point load vs distributed, what faces etc... It would make identifiing issues in your boundary conditions easier to spot.

matthiass March 20, 2012 15:46

Hi kmooney,

thanks for your reply.
The dimensions of the beam are (0.02m x 0.02m x 1m). Referring to the picture the bottom face (0.02m x 0.02m) is fixed (I called it clamp), on the left face (0.02m x 1m) a surface respectively distributed load of 10kN/m˛ is applied and the other faces are traction free.
I hope it gets clearer now.

Best regards
Matthias

matthiass April 17, 2012 06:49

Still no acceptable solution
 
5 Attachment(s)
Hello everybody,

until now I tried different things like changing the relative tolerance, large meshes, and so on. But I still get bizarre solutions. Because of that I decided to post again a bit more detailed.

I drew a sketch of my setup ->Balken.jpg.
I solved the problem analytically, with Comsol Multiphysics and with OpenFoam. For a beam with a length L of 4m I get similar results for all solution (comparison_beam_4x0.2x0.2.jpg, sigmaEq_beam_4x0.2x0.2.jpg). But for longer beams the results of OpenFoam are quite strange, e.g. for a beam with length L of 8m (comparison_beam_8x0.2x0.2.jpg, sigmaEq_beam_8x0.2x0.2.jpg) and I really don't know what I do wrong. The solution of OpenFoam for sigma shows a local peak inbetween the beam as illustrated in sigmaEq_beam_8x0.2x0.2.jpg, which leads to an inflection point in the displacement graph (see comparison_beam_8x0.2x0.2.jpg).

The only change between the 4m- and 8m-beam I made, was to change the length of the beam in blockMeshDict from 4 to 8 and the number of cells , so that the cells are cubes.

My course of action was:
blockMesh
decomposePar
mpirun -np 4 solidDisplacementFoam -parallel >& log
reconstructPar
Here is my code:
folder 0:
file D

Code:

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

dimensions      [0 1 0 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{
    left
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }

    right
    {
        type            tractionDisplacement;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value          uniform (0 0 0);
    }

    top
    {
        type            tractionDisplacement;
        traction        uniform ( 0 10000 0 );
        pressure        uniform 0;
        value          uniform (0 0 0);
    }

    bottom
    {
        type            tractionDisplacement;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value          uniform (0 0 0);
    }

    front
    {
        type            tractionDisplacement;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value          uniform (0 0 0);
    }

    back
    {
        type            tractionDisplacement;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value          uniform (0 0 0);
    }
}

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

folder constant
file mechanicalProperties

Code:

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

rho
{
    type        uniform;
    value      7854;
}

nu
{
    type        uniform;
    value      0.3;
}

E
{
    type        uniform;
    value      2e+11;
}

planeStress    no;


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

folder constant/polyMesh
file blockMeshDict

Code:

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

convertToMeters 1;

vertices
(
    (0 0 0)
    (4 0 0)
    (4 0.2 0)
    (0 0.2 0)
    (0 0 0.2)
    (4 0 0.2)
    (4 0.2 0.2)
    (0 0.2 0.2)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (120 6 6) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    left
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }
    right
    {
        type patch;
        faces
        (
            (2 6 5 1)
        );
    }
    top
    {
        type patch;
        faces
        (
            (3 7 6 2)
        );
    }
    bottom
    {
        type patch;
        faces
        (
            (1 5 4 0)
        );
    }
    front
    {
        type patch;
        faces
        (
            (0 3 2 1)
        );
    }
    back
    {
        type patch;
        faces
        (
            (4 5 6 7)
        );
    }
);

mergePatchPairs
(
);

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

folder system
file controlDict
Code:

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

application    solidDisplacementFoam;

startFrom      startTime;

startTime      0;

stopAt          endTime;

endTime        1;

deltaT          1;

writeControl    timeStep;

writeInterval  1;

purgeWrite      0;

writeFormat    ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision  6;

graphFormat    raw;

runTimeModifiable true;


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

file decomposeParDict
Code:

/*---------------------------------------------------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.3                                  |
|  \\  /    A nd          | Web:      http://www.openfoam.org              |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version        2.0;
    format          ascii;

    root            "";
    case            "";
    instance        "";
    local          "";

    class          dictionary;
    object          decomposeParDict;
}

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


numberOfSubdomains 4;

method          simple;

simpleCoeffs
{
    n              (4 1 1);
    delta          0.001;
}

hierarchicalCoeffs
{
    n              (1 1 1);
    delta          0.001;
    order          xyz;
}

metisCoeffs
{
    processorWeights
    (
        1
        1
        1
    );
}

manualCoeffs
{
    dataFile        "";
}

distributed    no;

roots
(
);


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

file fvSchemes
Code:

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

d2dt2Schemes
{
    default        steadyState;
}

ddtSchemes
{
    default        Euler;
}

gradSchemes
{
    default        leastSquares;
    grad(D)        leastSquares;
    grad(T)        leastSquares;
}

divSchemes
{
    default        none;
    div(sigmaD)    Gauss linear;
}

laplacianSchemes
{
    default        none;
    laplacian(DD,D) Gauss linear corrected;
    laplacian(DT,T) Gauss linear corrected;
}

interpolationSchemes
{
    default        linear;
}

snGradSchemes
{
    default        none;
}

fluxRequired
{
    default        no;
    D              yes;
    T              no;
}


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

file fvSolution
Code:

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

solvers
{
    "(D|T)"
    {
        solver          GAMG;
        tolerance      1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels    1;
    }
}

stressAnalysis
{
    compactNormalStress yes;
    nCorrectors    500000;
    D              1e-06;
}


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

I hope anyone can help me to find my mistake!
Thank you in advance
Matthias

lovecraft22 April 17, 2012 07:17

I don't think I could be of any help to you but probably you'll have more chances if you upload your latest case directory…

matthiass April 17, 2012 07:37

latest test case
 
1 Attachment(s)
Thank you for this advice!
Here is my test case, but unfortunately my solution is larger than 97KB.

Best regards,
Matthias

lovecraft22 April 17, 2012 09:33

Any specific reason why you changed your fvSolution file at the bottom if compared to the plateHole tutorial?

matthiass April 17, 2012 10:01

Yes, I did this because with the configuration of the plate hole tutorial my setup didn't converge due to too less iterations. I found this in the following thread http://www.cfd-online.com/Forums/ope...sanalysis.html
mentioned by bigphil in his 3rd post.
Quote:

Originally Posted by bigphil (Post 346473)
...
The nCorrectors should be much larger, I would typically set it to 5000. This will ensure that the model is converged for each time step.
The reason nCorrectors is set to 1 in the tutorial case is because the steady state results are achieved through multiple time steps (where the results are not converged for intermediate time-steps).

Philip


preibie April 17, 2012 10:37

1 Attachment(s)
Hi,

I take your case and find some strange things. You don't have any time loop you only use the gamg solver... I use your blockMeshDict and set a new case with the solidEquilibriumDisplacementFoam solver. In this are the stresses calculated to. I calculate an displacement of 0,03692 m.

I'am very interested in the comparison with the other FEM solver ant the analytic solution can you send me the picture with the new results ore the data files.

Stefan

bigphil April 17, 2012 11:24

Hi Matthiass,

I tried your case (8m beam) and I was able to get displacements much closer to the analytical. I think you just need to solve to a finer tolerance, I would change the D tolerance to 1e-7 or 1e-8 or maybe smaller.
i.e.
Code:

stressAnalysis
{
    compactNormalStress yes;
    nCorrectors    500000;
    D              1e-08;
}

This will probably take quite a few outer iterations but it should be better.

Philip

matthiass April 17, 2012 12:23

that's it
 
1 Attachment(s)
Hello,

thanks a lot for all those answers.

@ Stephan: I added the diagram with the solidEquilibriumDisplacementFoam and the result looks much better. I already tried the solidEquilibriumDisplacementFoam-solver but with more ncorrectors and smaller relative tolerance. But the result wasn't that good and so I tried to come along with solidDisplacementFoam-solver.

@ Philip: Did you only change the tolerance of D? Because I already tried this a view days ago and it took really long to solve. After a while the Initial and the Final residual showed the same solution (slightly under 1e-6) with No Iterations 0 for all three directions. This took again a while and the solution was not really better than before. At the moment I try it again with a tolerance of 1e-8. How long did the simulation take on your computer? Maybe my computer is too slow but he took over an hour. With Comsol it took a view seconds and with the setup of Stephan it took a view minutes.

In any case I have to get a better understanding of that solver!

Best regards,
Matthias

bigphil April 17, 2012 12:41

Quote:

Originally Posted by matthiass (Post 355199)
@ Philip: Did you only change the tolerance of D? Because I already tried this a view days ago and it took really long to solve. After a while the Initial and the Final residual showed the same solution (slightly under 1e-6) with No Iterations 0 for all three directions.

Hi Matthias,

You need to change the solver tolerance (GAMG or PCG or whatever you use) and also change the D tolerance. The solver tolerance should be tighter than the D tolerance otherwise the solver will do zero iterations as you have seen.
Use fvSolution looks like this:
Code:

    D
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance      1e-10;
        relTol          0.1;
    }
}

stressedFoam
{
    nCorrectors    200000;
    D              1e-09;
}

But yes this is very slow (~1 hour) and it needs many outer correctors.
This is because this is a bending dominated case and bending is strongly coupled in x-y-z directions. Since OpenFOAM solves the equations in a segregated manner (i.e. x then y then z and repeat) this will cause slow convergence in cases that are strongly coupled between directions.

The solidEquilibriumDisplacement solver may be a faster due to its use of an acceleration factor but I am not sure.

A block coupled solid solver would REALLY help in this case, maybe in time…


Philip

matthiass April 17, 2012 12:55

Hi Philip,

thank you very much! Now I got it! You helped me a lot ;-)

Best regards,
Matthias

sail April 17, 2012 17:39

I've seen that you have swithced of the plane stresses. (mechanichal proprietis file)


althought I agree that they might be really small in a case like this, maybe it can be worth a try to run a simulation bringing the plane stresses back on.

bigphil April 17, 2012 18:07

Quote:

Originally Posted by sail (Post 355244)
I've seen that you have swithced of the plane stresses. (mechanichal proprietis file)


althought I agree that they might be really small in a case like this, maybe it can be worth a try to run a simulation bringing the plane stresses back on.

Hmnn planeStress should only be set to yes for 2D thin cases.
planeStress should be set to no for all 3D cases,
as outlined in H. Jasak, H. G. Weller. Application of the finite volume method and unstructured meshes to linear elasticity. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 2000; 48:267–287.


Philip

preibie April 18, 2012 03:24

Hallo,

@bigphil & @matthiass

your calculate the beam with a high number of corrector steps and with a low number of "time" steps. In my opinion you should try to use only a small number of corrector steps may by 2 till 10. You have to know that OpenFoam is not an pure FEM solver. When you use the time step formulation it is possible to look for the results (paraview or shell). A small question: when do you think a calculation is converged? Not when the errors are below a limit like 10⁻⁶ or so. A problem is solved when the physic values are constant and errors are low! Did you look at the results after 2000 corrector steps?

Stefan

matthiass April 19, 2012 14:02

Hi Stefan,

thanks for your reply.
Until your post I didn't know that.

How can I formulate an abort criterion that stops when the physic values are constant and errors are low? Because in the case you sent me, all the time steps are calculated, even when the results were already well.
I think, I haven't already understood what happens if I use the time step formulation.

Best regards,
Matthias

preibie April 20, 2012 05:30

I don't know exactly how do you can create an abort criteria. I do it this way. I write me to the values I'm interested in an log file. Then I use foamLog to extract all values and plot the value in gnuplot. For all this you can make a small script.

Stefan

bigphil April 20, 2012 05:48

Stefan is right that the residual is not necessary a measure of the convergence of the physics (although it is normally a good measure).

Sometimes I will use a different residual based on when the displacement has stopped changing. If you create a variable before the momentum loop:
Code:

scalar relativeResidual = GREAT;
Then calculate this relative residual at the end of the momentum loop (after the solve function):
Code:

{
    scalarField magD = mag(D.internalField() - D.oldTime().internalField());

    forAll(magD, cellI)
    {
        if (magD[cellI] < SMALL)
        {
            magD[cellI] = SMALL;
        }
    }

    relativeResidual =
        gMax
        (
            mag
            (
                D.internalField()
              - D.prevIter().internalField()
            )
          /magD
        );
}

And then change the loop criteria:
Code:


          //solverPerf.initialResidual() > convergenceTolerance
          relativeResidual > convergenceTolerance

This is normally a good measure of the convergence of the physics. However, in the case of bending dominated models (as in this case), the segregated finite volume approach take a very long time to converge the physics (maybe hundreds of thousands of iterations).
Finite element solvers typically use a block coupled solution method so it has no problem with bending dominated cases (except it uses lots more RAM). The discretisation of the finite volume terms for a block coupled solver is more complicated but it is possible, and may be available in the coming years.

Philip

matthiass April 20, 2012 12:44

Hi,
thank you to both of you for all these answers and have a nice weekend!
Best regards,
Matthias


All times are GMT -4. The time now is 01:29.