[matthiass]

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:

/*--------------------------------*- 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]

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]

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]

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

/*--------------------------------*- 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

/*--------------------------------*- 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

    
/*--------------------------------*- 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

/*--------------------------------*- 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

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  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

/*--------------------------------*- 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

/*--------------------------------*- 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]

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]

Thank you for this advice!
Here is my test case, but unfortunately my solution is larger than 97KB.

Best regards,
Matthias

[lovecraft22]

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

[matthiass]

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.

Originally Posted by bigphil

...
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]

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]

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.

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

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

Philip

[matthiass]

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]

Originally Posted by matthiass

@ 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:

    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]

Hi Philip,

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

Best regards,
Matthias

[sail]

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]

Originally Posted by sail

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]

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]

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]

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]

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:

scalar relativeResidual = GREAT;

Then calculate this relative residual at the end of the momentum loop (after the solve function):

{
    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:

 
           //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]

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