
[Sponsors] 
[stressAnalysis] Problems with simple test case 

LinkBack  Thread Tools  Search this Thread  Display Modes 
March 20, 2012, 09:01 
Problems with simple test case

#1 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 14 
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, 09:48 

#2 
Senior Member
Kyle Mooney
Join Date: Jul 2009
Location: San Francisco, CA USA
Posts: 323
Rep Power: 17 
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, 14:46 

#3 
New Member
Join Date: Mar 2012
Posts: 23
Rep Power: 14 
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: 14 
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: 14 
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: 14 
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.cfdonline.com/Forums/ope...sanalysis.html
mentioned by bigphil in his 3rd post. Quote:


April 17, 2012, 10:37 

#9 
Member
Stefan
Join Date: Jun 2009
Posts: 67
Rep Power: 17 
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 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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: 14 
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  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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: 14 
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: 16 
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  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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: 17 
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: 14 
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: 17 
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 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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: 14 
Hi,
thank you to both of you for all these answers and have a nice weekend! Best regards, Matthias 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Zetaf turbulence model test case  MasAmi  OpenFOAM Running, Solving & CFD  1  March 3, 2014 18:31 
T3c2 flat plate test case  rajeshamech  Main CFD Forum  0  November 21, 2013 05:21 
Help Required With Simple Test Case  steph79  OpenFOAM PreProcessing  4  August 3, 2010 07:45 
Interfoam Droplet under shear test case  adona058  OpenFOAM Running, Solving & CFD  3  May 3, 2010 18:46 
Combustion Test Case  A.S.  Main CFD Forum  1  May 31, 2005 09:22 