Question on how strain and stress is calculated in TetrahedronFEMForceField

[nhnhan92]

Hi everyone,

Here I would like to raise a question on the code of TetrahedronFEMForceField, especially for the algorithm to calculate stress and strain. I have been going through the code (computeVonMisesStress()) and as far as I understand, the calculation method is as follows:

1. Calculate the displacement gradient = shapefunction * Displacement
2. Calculate strain using Green Strain Tensor (either full or short form): epsilon = 1/2(displacement gradient + displacement gradient^T)
3. Calculate stress = (equal to the sum of diagonal elements of strain)

From the above procedure, I have two following concerns that I hope someone can clarify them in more detail:
1.

Mat33 strain = Real(0.5)*(gradU + gradU.transposed());`
 for (Index i = 0; i < 3; i++)
vStrain[i] = strain[i][i];
vStrain[3] = strain[1][2];
vStrain[4] = strain[0][2];
vStrain[5] = strain[0][1];

The vector vStrain contains strain components. The first three components (i = 1,2,3) are nothing to complain about. Regarding the last three components, is there any special reason to assign them like this? According to the strain_displcement_matrix, I think they are supposed to be like this:

vStrain[3] = strain[0][1];   \\\ tangential strain_xy
vStrain[4] = strain[1][2];   \\\ tangential strain_yz
vStrain[5] = strain[0][2];   \\\ tangential strain_xz

2. Strain from step 2 can also be calculated using the following equation: epsilon = strain_displcement_matrix * Displacement. However, I tried this method, but its strain outcome differs from that of the above methods. The weird thing is the calculated strain of the first element (elementIndex = 0) from both ways is equal. I am trying to figure out what is wrong. I guess something is wrong with the "Shapefunction - elemShapeFun[el]."

I hope someone might shed light on these problems. Thank you for your help in advance.

[hugtalbot]

Hey @nhnhan92

Thank you for your question, good to see you back in the code 🔥

The class TetrahedronFEMForceField implements a linear elastic model for small deformations but possibly large displacements. In such "small deformation" case, the Green Lagrange strain tensor can be linearized and noted ϵ:

Even if this model is only valid for small deformation, a co-rotational approach may be used to extract the rotation from the strain computation (see the option data method=large or method=polar).

Now about your questions:

1. the fact is that all tensors (matrices) are stored using the Voigt notation which reducing the size of these tensors due to their symmetric properties. But your understanding is correct.
2. yes it is correct :

I am not sure to understand exactly the problem: from what I understand you are surprised by the fact that the _**per element**_ strain is not equal to the previously computed strain tensor. Is this correct? I think this is due to the fact that in the physics all we need for the computation is the symmetric part of ϵ which differs from ϵ. Could you confirm this?

FYI : about the draw of the VonMises stress, a PR is currently open to update it

Note that to get more into deep on these topics, we now have a new training course on continuum mechanics and the Finite-Element Method (FEM)

Hope this helps.

[nhnhan92]

Dear @hugtalbot,

Thank you for your reply. I think I did not make a clear statement on point No. 2. In the code of TetrahedronFEMForceField, the strain tensor (i.e., variable "vStrain" in the code below) for each element is calculated as shown below:

Mat33 strain = Real(0.5)*(gradU + gradU.transposed());
 for (Index i = 0; i < 3; i++)
       vStrain[i] = strain[i][i];
     vStrain[3] = strain[1][2];
     vStrain[4] = strain[0][2];
      vStrain[5] = strain[0][1];`

I tried to compare the above strain output with those computed using the equation: epsilon = strain_displcement_matrix * Displacement. The below code shows how I execute this equation:

template<class DataTypes>
inline void TetrahedronFEMForceField<DataTypes>::computeStrain()
{
    if(this->d_componentState.getValue() == ComponentState::Invalid)
        return ;
    typename core::behavior::MechanicalState<DataTypes>* mechanicalObject;
    this->getContext()->get(mechanicalObject);
    const VecCoord& X = mechanicalObject->read(core::ConstVecCoordId::position())->getValue();
    const VecCoord X0=this->mstate->read(core::ConstVecCoordId::restPosition())->getValue();

    helper::WriteAccessor<Data<type::vector<type::Vec6f> > > vstrain_total =  _totalstrain;

    typename VecElement::const_iterator it;
    Element index = *it;
    Index elementIndex = el;
    for(it = _indexedElements->begin(), el = 0 ; it != _indexedElements->end() ; ++it, ++el)
    {
                Index a = (*it)[0];
                Index b = (*it)[1];
                Index c = (*it)[2];
                Index d = (*it)[3];
                // Rotation matrix (deformed and displaced Tetrahedron/world)
                Transformation R_0_2;
                Displacement D;
                computeRotationLarge( R_0_2, X, a,b,c);

                type::fixed_array<Coord,4> rotatedInitialElements;
                rotatedInitialElements[0] = R_0_2*(X)[a];
                rotatedInitialElements[1] = R_0_2*(X)[b];
                rotatedInitialElements[2] = R_0_2*(X)[c];
                rotatedInitialElements[3] = R_0_2*(X)[d];
                // positions of the deformed and displaced Tetrahedron in its frame
                type::fixed_array<Coord,4> deforme;
                for(int i=0; i<4; ++i)
                    deforme[i] = R_0_2*X[index[i]];

                deforme[1][0] -= deforme[0][0];
                deforme[2][0] -= deforme[0][0];
                deforme[2][1] -= deforme[0][1];
                deforme[3] -= deforme[0]; 
                // displacement
                D[0] = 0;
                D[1] = 0;
                D[2] = 0;
                D[3] = _rotatedInitialElements[elementIndex][1][0] - deforme[1][0];
                D[4] = 0;
                D[5] = 0;
                D[6] = _rotatedInitialElements[elementIndex][2][0] - deforme[2][0];
                D[7] = _rotatedInitialElements[elementIndex][2][1] - deforme[2][1];
                D[8] = 0;
                D[9] = _rotatedInitialElements[elementIndex][3][0] - deforme[3][0];
                D[10] = _rotatedInitialElements[elementIndex][3][1] - deforme[3][1];
                D[11] =_rotatedInitialElements[elementIndex][3][2] - deforme[3][2];
                StrainDisplacement	J;
                computeStrainDisplacement(J, rotatedInitialElements[0], rotatedInitialElements[1], rotatedInitialElements[2], rotatedInitialElements[3]);
                if(_updateStiffnessMatrix.getValue())
                {
                    strainDisplacements[elementIndex][0][0]   = ( - deforme[2][1]*deforme[3][2] );
                    strainDisplacements[elementIndex][1][1] = ( deforme[2][0]*deforme[3][2] - deforme[1][0]*deforme[3][2] );
                    strainDisplacements[elementIndex][2][2]   = ( deforme[2][1]*deforme[3][0] - deforme[2][0]*deforme[3][1] + deforme[1][0]*deforme[3][1] - deforme[1][0]*deforme[2][1] );

                    strainDisplacements[elementIndex][3][0]   = ( deforme[2][1]*deforme[3][2] );
                    strainDisplacements[elementIndex][4][1]  = ( - deforme[2][0]*deforme[3][2] );
                    strainDisplacements[elementIndex][5][2]   = ( - deforme[2][1]*deforme[3][0] + deforme[2][0]*deforme[3][1] );

                    strainDisplacements[elementIndex][7][1]  = ( deforme[1][0]*deforme[3][2] );
                    strainDisplacements[elementIndex][8][2]   = ( - deforme[1][0]*deforme[3][1] );

                    strainDisplacements[elementIndex][11][2] = ( deforme[1][0]*deforme[2][1] );
                }
                type::VecNoInit<6,Real> JtD;
                JtD = J.multTranspose(D);

                Coord A = (X0)[b] - (X0)[a];
                Coord B = (X0)[c] - (X0)[a];
                Coord C = (X0)[d] - (X0)[a];
                Coord AB = cross(A, B);
                Real volumes6 = fabs( dot( AB, C ) );

                JtD  /= (volumes6);
                vstrain_total[es] = JtD;
            } 

But they are not equal to each other. I hope someone can point out any mistakes or misunderstandings here.
Thank you

[hugtalbot]

Hi @nhnhan92

Getting back to this topic, it seems to me that the two approaches are computing two different strains indeed:

• yours uses the full (global) displacement information U
• while the TetrahedronFEMForceField code extract rotations to compute the new displacement "D[]" in a co-rotational way
  with two different measurements of the displacement, you get two different measurements of strain

Does this make sense?
We could happily rediscuss directly if you want

[nhnhan92]

Dear @hugtalbot ,

Sound reasonable, I will give your suggestion a try and update the result later.

[hugtalbot]

Great keep us posted 👍

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@nhnhan92
Feedback has been given to you by the project reviewers, however we have not received a response from you. Without further news in the coming weeks, this discussion will be automatically closed in order to keep this forum clean and fresh 🌱 Thank you for your understanding