[Sargam05]

Hello everyone,

I am doing compression test in two-phase microstrucutre geometry to calculate stresses and strains in different loading conditions. In the compression test, I have assumed uni-axial 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 stress-strain 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 two-phase (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 non-percolating (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 two-phase 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

[bigphil]

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 OpenFOAM-1.6-ext.

Hope it helps,
Philip

[Sargam05]

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 (OpenFOAM-ext) 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 stress-average 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.6511E-11

When I am solving for solid cube geometry, I got sigmaAv=10, and epsilonyy=4.6511E-11 (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.6511E-11 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

[bigphil]

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

[Sargam05]

Hi Philip,

Thank you for the quick reply! As you suggested, now I am running program in both OpenFOAM-ext with fixedDisplacementZeroShear BC and OpenFOAM-1.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 two-phase microstrucutre as attachment.

Best regards,
Sangeeta

[Sargam05]

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

[bigphil]

Sargam,

Originally Posted by Sargam05

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_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);
}
}

Your directionMixed refValue seems way too large, these are the normal displacements you are enforcing.
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

[Sargam05]

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.6511E-11 0);
refGradient uniform (0 0 0);
valueFraction uniform (0 0 0 1 0 0);
}
electron_maxY
{
type directionMixed;
refValue uniform (0 -4.6511E-11 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 y-direction is 4.6511E-11, 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 y-direction. 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.6511E-11 0);
}
electron_maxY
{
type fixedValue;
value uniform (0 -4.6511E-11 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

[Sargam05]

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

[bigphil]

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

                //DEqn.setComponentReference(1, 0, vector::X, 0);                                                                                             
                //DEqn.setComponentReference(1, 0, vector::Z, 0);

Also newStressedFoam allows you to specify the patch and face to fix (see plateHole/system/fvSolution):

    componentReference
    (
        {
            patch down;
            face  0;
            direction z;
            value 0;
        }
    );

By the way, there are a number of people in my research group who perform similar micro-structural simulations on representative volume elements and they normally define the effective E using a homogenisation technique:
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

[Sargam05]

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

[Sargam05]

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

[bigphil]

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

[Sargam05]

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 two-phase 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

[bigphil]

Originally Posted by Sargam05

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 think that in this case the effective stress essentially means the average stress, and using an approach similar to homogenisation is probably the correct method. For example, the stresses on the top patch may not be a good representation of the stress in the entire RVE so using homogenisation you can get a good average.

Originally Posted by Sargam05

I have refine my two-phase 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.

Yep I would increase the number of iterations - I have often had simulations requiring over 100,000 iterations or even over 1 million.

Originally Posted by Sargam05

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.

There should be no shear stresses on the boundary patches if you are only specifying normal stresses - if your model is not converged then it may show erroneous shear stresses.
But shear stresses in the interior of the model are to be expected given the complex geometry.

Originally Posted by Sargam05

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.

As mentioned before, when you use fixedValue this fixes the block from contracting/expanding in the tangential direction so this will generate boundary shear stresses and produce a stiffer model (there will be plane strain conditions in the centre of the model).

So I would work on trying to get the traction model to converge or try get the directionMixed model to converge.

Philip

[Sargam05]

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 strain-stress 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

[bigphil]

Hi Sangeeta,

Apologies for the delay in replying.

Originally Posted by Sargam05

I have also tried following directionMixed boundary condition for a solid cube geometry to check validity for the boundary conditions:
electron_minY
{
type directionMixed;
refValue uniform (0 0 0);
refGradient uniform (0 10 0);
valueFraction uniform (0 0 0 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?

With the directionMixed BC above you are forcing a surface normal gradient of 10! - I really don't think this is what you intended - this is essentially straining the boundary to 1000%!

I would recommend re-reading 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