
[Sponsors] 
March 20, 2012, 10:01 
Problems with simple test case

#1 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 FEsoftware 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 Dfile, 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); } } // ************************************************************************* // 

March 20, 2012, 10:48 

#2 
Senior Member
Kyle Mooney
Join Date: Jul 2009
Location: Amherst, MA USA  San Diego, CA USA
Posts: 268
Rep Power: 8 
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.


March 20, 2012, 15:46 

#3 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 

April 17, 2012, 06:49 
Still no acceptable solution

#4 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 8mbeam 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: blockMeshHere 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); } } // ************************************************************************* // 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; // ************************************************************************* // 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 ( ); // ************************************************************************* // 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; // ************************************************************************* // 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 ( ); // ************************************************************************* // 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; } // ************************************************************************* // 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 { "(DT)" { solver GAMG; tolerance 1e06; relTol 0.01; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 20; agglomerator faceAreaPair; mergeLevels 1; } } stressAnalysis { compactNormalStress yes; nCorrectors 500000; D 1e06; } // ************************************************************************* // Thank you in advance Matthias 

April 17, 2012, 07:17 

#5 
Senior Member

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…


April 17, 2012, 07:37 
latest test case

#6 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
Thank you for this advice!
Here is my test case, but unfortunately my solution is larger than 97KB. Best regards, Matthias 

April 17, 2012, 09:33 

#7 
Senior Member

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


April 17, 2012, 10:01 

#8  
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 Someone uses the OpenFOAM for StressAnalysis?
mentioned by bigphil in his 3rd post. Quote:


April 17, 2012, 10:37 

#9 
Member
Stefan
Join Date: Jun 2009
Posts: 67
Rep Power: 8 
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 

April 17, 2012, 11:24 

#10 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 565
Rep Power: 19 
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 1e7 or 1e8 or maybe smaller. i.e. Code:
stressAnalysis { compactNormalStress yes; nCorrectors 500000; D 1e08; } Philip 

April 17, 2012, 12:23 
that's it

#11 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 solidEquilibriumDisplacementFoamsolver but with more ncorrectors and smaller relative tolerance. But the result wasn't that good and so I tried to come along with solidDisplacementFoamsolver. @ 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 1e6) 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 1e8. 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 

April 17, 2012, 12:41 

#12  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 565
Rep Power: 19 
Quote:
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 1e10; relTol 0.1; } } stressedFoam { nCorrectors 200000; D 1e09; } This is because this is a bending dominated case and bending is strongly coupled in xyz 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 

April 17, 2012, 12:55 

#13 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
Hi Philip,
thank you very much! Now I got it! You helped me a lot ;) Best regards, Matthias 

April 17, 2012, 17:39 

#14 
Senior Member
Vieri Abolaffio
Join Date: Jul 2010
Location: Always on the move.
Posts: 308
Rep Power: 7 
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. 

April 17, 2012, 18:07 

#15  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 565
Rep Power: 19 
Quote:
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 

April 18, 2012, 03:24 

#16 
Member
Stefan
Join Date: Jun 2009
Posts: 67
Rep Power: 8 
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 

April 19, 2012, 14:02 

#17 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
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 

April 20, 2012, 05:30 

#18 
Member
Stefan
Join Date: Jun 2009
Posts: 67
Rep Power: 8 
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 

April 20, 2012, 05:48 

#19 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 565
Rep Power: 19 
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; 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 ); } Code:
//solverPerf.initialResidual() > convergenceTolerance relativeResidual > convergenceTolerance 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 

April 20, 2012, 12:44 

#20 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 5 
Hi,
thank you to both of you for all these answers and have a nice weekend! Best regards, Matthias 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Interfoam Droplet under shear test case  adona058  OpenFOAM Running, Solving & CFD  3  May 3, 2010 18:46 
Simple problems with ansys CFX  Luk_Fiz  CFX  11  October 26, 2009 18:17 
3D laminar flow test case  Praveen. C  Main CFD Forum  2  March 6, 2008 05:28 
Flat Plate Transition Test Case  Adam  CFX  1  June 23, 2005 23:12 
Durham test case  SAM  FLUENT  0  August 16, 2004 05:01 