
[Sponsors] 
OpenFOAM solution is diverging for stress analysis in twopahse microstructure. 

LinkBack  Thread Tools  Display Modes 
April 6, 2013, 00:30 
OpenFOAM solution is diverging for stress analysis in twopahse microstructure.

#1 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hello everyone,
I am doing compression test in twophase microstrucutre geometry to calculate stresses and strains in different loading conditions. In the compression test, I have assumed uniaxial loading on top (maxY) and bottom (minY) surfaces and other surfaces have zero traction (i.e. free surfaces). I have first validated OpenFOAM solution with analytical solution. I have used following displacement boundary conditions for the stressstrain analysis in solid cube geometry: /** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.6   \\ / A nd  Web: http://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 { minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } minY { type tractionDisplacement; traction uniform ( 0 10 0 ); pressure uniform 0; value uniform (0 0 0); } maxY { type tractionDisplacement; traction uniform ( 0 10 0 ); pressure uniform 0; value uniform (0 0 0); } minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } OpenFOAM solution for the above displacement boundary conditions, matches with analytical solution. I have generated a twophase (spherical particle  pore) microstructure. I named spherical particles as electrons. I am not considering any pore or air pressure for the analysis. Generated microstucture geometry have some percolating pore (pore) and nonpercolating (npore), therefore I have two more patches here. I have used same boundary conditions which I have validated for the solid cube geometry (which had same size as microstructure). The boundary conditions for twophase geometry are: /** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.7.1   \\ / A nd  Web: www.OpenFOAM.com   \\/ 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 { electron_minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minY { type tractionDisplacement; traction uniform ( 0 10 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxY { type tractionDisplacement; traction uniform ( 0 10 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_pore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_npore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } } But I am getting diverge solutions with "nan" values. i do not know why this is happening. Could anyone have any idea why this problem is happening? What boundary conditions should I use for the compression test so that I could have converged solution? I appreciate if anyone can give me some inputs/suggestions. Thanks in advance. Best regards, Sangeeta 

April 6, 2013, 14:38 

#2 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Hi Sangeeta,
I presume the models are static i.e. you are ignoring inertia. If this is the case, then your model setup is very unstable as you are only using Neumann (traction) boundary conditions, you should have at least one Dirichlet (fixed displacement) condition. The simple model probably achieved convergence because it was simple  one material, probably simple orthogonal mesh. So maybe try change the top and bottom boundaries to fixedValue with displacement equivalent to that induced by the tractions. Secondly, are you assigning different mechanical properties to each phase? Are the stiffnesses very different? If the stiffnesses are very different then the tractions at the interfaces between the phases will be poorly approximated using the standard OpenFOAM discretisation. A procedure like that proposed by Tukovic et al is required as implemented in the solidInterface class in the solidMechanics branch of OpenFOAM1.6ext. Hope it helps, Philip Last edited by bigphil; April 6, 2013 at 14:38. Reason: typo 

April 6, 2013, 15:50 

#3 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Thank you so much for the reply! I am only considering mechanical properties of solid phase (spherical particles). This is a linear elastic steady state problem. I have tried SolidDispalcementFoam (OpenFOAM 1.7.1) and elasticSolidFoam (OpenFOAMext) to validate the case for solid cube geometry and the solution given by these solvers matches with analytical solution. In this problem, I need to calculated average stressaverage strain relationship and effective Young's modulus. The average stresses and stains need to measure for random positions of electrons in the domain, therefore I cannot assume one or two symmetry planes (because I want all electrons in random places, i.e. no symmetry). In addition, I need average values of stresses and strains which cannot be accurately measured without defining the all surfaces (i.e. without symmetry planes). I am not solving anything for pore phase (i.e. not considering any pore pressure or species). Following mechanical properties have used for the validation: FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object mechanicalProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // rho rho [ 1 3 0 0 0 0 0 ] 5163; nu nu [ 0 0 0 0 0 0 0 ] 0.319; E E [ 1 1 2 0 0 0 0 ] 215e+09; planeStress no; As you suggested to change the top and bottom boundaries to fixedValue with displacement equivalent to that induced by the tractions. Do you mean that I should analytically calculate strain and then change in length (which is displacement), and use this value for fixed displacement BC. Here is the analytical calculation based on the mechanical properties. stress = E*strain load/area = E*change in length/original length 10/(1X1) = (215 E+09)* change in length/1 change in length = 4.6511E11 When I am solving for solid cube geometry, I got sigmaAv=10, and epsilonyy=4.6511E11 (which is change in length) Do you mean that I should used the change in length value for fixedValue displacement boundary condition on top and bottom surfaces of microstructure as: type fixedValue; value uniform (0 4.6511E11 0); and other surfaces would be free traction surfaces (free surfaces) including electron_to_pore and electron_to_npore. Please let me know what did you mean? Do I need to use traction free BC for electron_to_pore and electron_to_npore? Best regards, Sangeeta 

April 6, 2013, 16:11 

#4 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Hi Sangeeta,
essentially yes, apply fixedValue to the top and bottom with the values you have calculated. But remember that when you use fixedValue you may get significant shear stresses on these boundaries, so you could use a BC such as fixedDisplacementZeroShear which applies a displacement in the normal direction and enforces zero shear stress in the tangential direction (equivalent to directionMixed). If you do use fixedValue and you do have significant shear stresses then your model is no longer uniaxial so your analytical solution does not apply i.e. stress = YoungsModulus*strain is only strictly valid for uniaxial test cases. Also, I am not really sure I understand what your geometry looks like, could you post an image/schematic? Philip 

April 6, 2013, 20:38 

#5 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Thank you for the quick reply! As you suggested, now I am running program in both OpenFOAMext with fixedDisplacementZeroShear BC and OpenFOAM1.7.1 with directionMixed BC on top and bottom surfaces respectively. As I am applying uniaxial or normal load on minY and maxY directions, I have considered following BDc in OpenFOAM1.7.1 FoamFile { version 2.0; format ascii; class volVectorField; object D; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 0 0 0 0 0]; internalField uniform (0 0 0); boundaryField { electron_minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minY { type directionMixed; refValue uniform (0 10 0); refGradient uniform (0 0 0); valueFraction uniform (0 0 0 1 0 0); } electron_maxY { type directionMixed; refValue uniform (0 10 0); refGradient uniform (0 0 0); valueFraction uniform (0 0 0 1 0 0); } electron_minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_pore { type tractionDisplacement; traction uniform ( 0.00001 0.00001 0.00001 ); pressure uniform 0; value uniform (0 0 0); } electron_to_npore { type tractionDisplacement; traction uniform ( 0.00001 0.00001 0.00001 ); pressure uniform 0; value uniform (0 0 0); } } Do you think I am using directionMixed BC is fine for the case? Also please find the twophase microstrucutre as attachment. Best regards, Sangeeta 

April 7, 2013, 00:09 

#6 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
I am still getting diverge solutions even after using above boundary condition. Do you have any idea what I am doing wrong here? Best regards, Sangeeta 

April 7, 2013, 08:52 

#7  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Sargam,
Quote:
Also, your values for traction are very small and are essentially zero. Have you got your displacement and traction values mixed up? By the way, to begin with I would try the case with fixedValue BCs instead of directionMixed as fixedValue are more stable. Could you point out where each of your boundaries are on the model geometry? as it is not clear where each boundary is. Philip 

April 7, 2013, 12:15 

#8 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Thanks for the reply. Mistakenly I have sent wrong 0/D file to you. the original boundary conditions I am using are: FoamFile { version 2.0; format ascii; class volVectorField; object D; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 0 0 0 0 0]; internalField uniform (0 0 0); boundaryField { electron_minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minY { type directionMixed; refValue uniform (0 4.6511E11 0); refGradient uniform (0 0 0); valueFraction uniform (0 0 0 1 0 0); } electron_maxY { type directionMixed; refValue uniform (0 4.6511E11 0); refGradient uniform (0 0 0); valueFraction uniform (0 0 0 1 0 0); } electron_minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_pore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_npore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } } Here you could see that I am uniaxial displacement in ydirection is 4.6511E11, and valueFraction is (0 0 0 1 0 0) i.e. sigmayy=1. I have also attached the detailed geometry with boundaries in it. The geometry shows that electron_to_pore boundary in outside boundaries between pore and electron. electron_to_npore refers to non percolating pore (npore) space inside of the geometry. I considered electron_to_pore and electron_to_npore as free surfaces. The load is coming only in normal ydirection. I am getting diverge solution when i am using these BCs. As you suggested, I tried a fixedValue case first with following BCs: FoamFile { version 2.0; format ascii; class volVectorField; object D; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 0 0 0 0 0]; internalField uniform (0 0 0); boundaryField { electron_minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minY { type fixedValue; value uniform (0 4.6511E11 0); } electron_maxY { type fixedValue; value uniform (0 4.6511E11 0); } electron_minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_pore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_to_npore { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } } The solution I am getting is not matching to physical solution that mean the Young's modulus (E) I got from the solution is higher than the Young's modulus i have taken for the calculation. Theoretically, E should be decreased with increasing pore space in the geometry. Do you think it could be because of shear stresses which are generating due to fixedValue BC? Am I doing something wrong here? Please let me know what do you think? Best regards, Sangeeta 

April 7, 2013, 12:20 

#9 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Sorry Philip, I am not able to send the geometry with BCs here because it has
exceeded the desired size of the attachment. I am sending the geometry in your mail address. Thanks, Sangeeta 

April 7, 2013, 13:04 

#10 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Sangeeta,
OK thanks I understand your case better now. So the case works with fixedValue BCs but the apparent E is too high  this is to be expected because of the shears  the model will seem stiffer because you are constraining it at the top and bottom from expanding. Therefore you need directionMixed style BCs but the problem with using directionMixed for the top and bottom boundaries in your model is that the model is not constrained in the X and Z directions. For example if your model starts to move in the X/Z direction there is nothing to stop it and it will fly off to infinity and never converge. So you need to constrain the motion of the model in X and Z directions  maybe symmetryPlanes on mixX and minZ? Another option is to set a reference component, this essentially means that during the solution procedure we keep fixing the X/Z component of displacement to zero at a chosen cell/face so the model doesn't fly away. Have a look at solidDisplacementFoam (line 94/95): Code:
//DEqn.setComponentReference(1, 0, vector::X, 0); //DEqn.setComponentReference(1, 0, vector::Z, 0); Code:
componentReference ( { patch down; face 0; direction z; value 0; } ); E_eff = sigma_eff / epsilon_eff where the effective stress sigma_eff = (1/V) integral_over_volume sigma dV and the effective strain epsilon_eff = (1/V) integral_over_volume epsilon dV where the volume integrals are performed over the entire cube (not just the particles). Best regards, Philip 

April 8, 2013, 00:03 

#11 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Thanks a lot for the valuable inputs! I have compiled solidDisplacmentFoam solver with DEqn.setComponentReference(1, 0, vector::X, 0); DEqn.setComponentReference(1, 0, vector::Z, 0); and set fixed displacement BCs on top and bottom surfaces. Program is still running. This time solution is not giving infinite values but it seems that still calculated values different than analytical solution. You have mentioned about homogenisation technique for solving effective properties. Did you mean total mesh volume? Could you so kind to elaborate the technique? Best regards, Sangeeta 

April 10, 2013, 11:57 

#12 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
I have tried different methods you have suggested, and now I got it! Thank you so very much for the inputs and suggestions. Best regards, Sangeeta 

April 11, 2013, 05:38 

#13 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Hi Sangeeta,
Sorry for the delay in replying. That's good to hear it is working for you, how exactly did you solve your problem for future reference? For the homogenisation, the stress and strain are integrated over the entire cube volume (so not necessarily the mesh volume). Integrating the stress is straightforward as the stress is zero in the porous regions, but the strain in the porous regions is not zero and must be included in the integration. Philip 

April 11, 2013, 11:30 

#14 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Actually, I need to calculate average stresses and strains not effective stresses and strains for my geometry because I need to generate some data of material properties based on geometry. I have refine my twophase mesh geometry and now solution is not diverging for previous BCs (traction forces on top and bottom and other surfaces are traction free), but still I did not get converge solution yet. I set 1000 ncorrector for 1 iteration but solution did not converged. I think I need to run it for more number of iterations. But again the solution I got is showing higher values of E than actual used value even after setting reference values of x and z to zero in DEqn. It seems shear stresses and strains are still generating. For a solid cube geometry, when I used traction BCs (i.e. traction forces on top and bottom and other surfaces are traction free), the calculated average stress value was equal to the traction force value which I used for the analysis. I used fixed displacement conditions (fixed displacement has calculated analytically for some force) on normal to the top and bottom surfaces for a solid cube geometry. The calculated average stress value is not equal to the traction force which I have taken to calculate displacement. Theoretically, calculated stress value will be equal to considered value of force for both fixed displacement and traction force BCs. I do not know why this is happening. I appreciate if you have any suggestion for the problem. Thank you. Best regards, Sangeeta 

April 14, 2013, 15:18 

#15  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Quote:
Quote:
Quote:
But shear stresses in the interior of the model are to be expected given the complex geometry. Quote:
So I would work on trying to get the traction model to converge or try get the directionMixed model to converge. Philip 

April 15, 2013, 12:48 

#16 
Member
Sangeeta
Join Date: Jul 2012
Location: Kingston, Canada
Posts: 70
Rep Power: 6 
Hi Philip,
Thank you for the reply! I put 100 necorrectors and 800 iterations for the simulation but my model is running for 3 days. If I would set more number of iteration let 1 million, it will take very long to complete simulation. Should I reduce number of ncorrectors? I have not change equation in solidDisplacement foam. I have only added strain equation in .C file and strainstress tensors in calculateStress file. I do not know what is wrong with the program why it is taking too ling to run. I have also tried following directionMixed boundary condition for a solid cube geometry to check validity for the boundary conditions: dimensions [0 1 0 0 0 0 0]; internalField uniform (0 0 0); boundaryField { electron_minX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxX { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_minY { type directionMixed; refValue uniform (0 0 0); refGradient uniform (0 10 0); valueFraction uniform (0 0 0 0 0 0); } electron_maxY { type directionMixed; refValue uniform (0 0 0); refGradient uniform (0 10 0); valueFraction uniform (0 0 0 0 0 0); } electron_minZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } electron_maxZ { type tractionDisplacement; traction uniform ( 0 0 0 ); pressure uniform 0; value uniform (0 0 0); } } But solution did not match with the analytical results by using above BCs. I have read various openFOAM Forum threads about directionMixed BC but still I am little confused to use it. I set reference components to be zero gradient dUy is 10 (similar to traction boundary condition where I set traction in y direction as 10). What is wrong the BCs I have used? What do you think? Best regards, Sangeeta 

April 30, 2013, 16:18 

#17  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 620
Rep Power: 22 
Hi Sangeeta,
Apologies for the delay in replying. Quote:
I would recommend rereading the other posts on the forum describing directionMixed and trying a few simple cases to get an understanding of it. In brief, refValue is a displacement, refGradient is a gradient and the valueFraction specifies how these are "mixed". A valueFraction of (0 0 0 0 0 0) is the same as a fixedGradient BC. A valueFraction of (1 0 0 1 0 1) is the same as fixedValue. Note that a gradient is not the same as a traction, a gradient is more similar to a strain, but you can calculate the gradient based on a prescribed traction. Philip 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Adventure of fisrst openfoam installation on Ubuntu 710  jussi  OpenFOAM Installation  0  April 24, 2008 14:25 
IcoFoam parallel woes  msrinath80  OpenFOAM Running, Solving & CFD  9  July 22, 2007 02:58 
HELP! Diverging Solution for Incompressible Prob.  Bob  FLUENT  3  April 16, 2006 19:33 
Could anybody help me see this error and give help  liugx212  OpenFOAM Running, Solving & CFD  3  January 4, 2006 19:07 
help: analysis solution of gaswater sod problem  gooo  Main CFD Forum  0  April 19, 2005 01:05 