Principal stress components

[Raúl Aparicio Yuste]

Hello everyone,

It is the first time I am using the library Deal II. I am wondering, is there any way/example to get the principal stress components as an output of the simulation? Thank you very much in advance

Raul

[Jan Philipp Thiele]

If I remember my mechanics courses correctly they can be calculated in a post-processing step after obtaining stresses and deformations
by solving the corresponding problem, right?
If so, there are multiple ways.
The maybe not so nice but quick way would be to just calculate the quantities in VisIt or ParaView when loading the simulation results in there.

The nicer way is to use the DataPostprocessor in deal.II itself to calculate quantities depending on the solution when doing the output.

This is exemplarily used in

https://dealii.org/developer/doxygen/deal.II/step_32.html

And there was also a code (possibly in the gallery) that uses the DataPostprocessor to calculate Vorticity based on the solution to Navier-Stokes equations.

I hope this helps.

Philipp

[Raúl Aparicio Yuste]

Dear Philipp,

Thank you for your answer. As you mention, one can compute the principal stress components in different ways, for example, by computing the eigen stress values or by using the Mohr’s circle. I thought that maybe someone had already faced that problem before with Dealii, or that Dealii itself could give you that variable as an output. That would be great. 

In the meantime, I am going to look at the examples you suggest, thank you very much.

Raul

[Wolfgang Bangerth]

On 11/18/22 05:42, Raúl Aparicio Yuste wrote:
> Thank you for your answer. As you mention, one can compute the principal
> stress components in different ways, for example, by computing the eigen
> stress values or by using the Mohr’s circle. I thought that maybe someone had
> already faced that problem before with Dealii, or that Dealii itself could
> give you that variable as an output. That would be great.

deal.II doesn't already have this built in, but it would certainly be very
nice if you contributed whatever you come up with for inclusion!

(I will caution that one of the reasons deal.II doesn't have this is because
the stress depends on what formulation of elasticity you use. It is relatively
straightforward if you have linear isotropic elasticity because then you can
compute it from the strain and just the two Lame parameters. But if you have
some more complex, possibly nonlinear material, then the computation of the
stress tensor can be substantially more complicated -- and that makes it less
appealing to include in a generic library. But it would make for a nice
example in the documentation of DataPostprocessorVector, for example, in the
same way as the computation of the stress tensor is currently in the
documentation of DataPostprocessTensor.)

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[richard....@gmail.com]

Hi,

I just recently did something similar and as has been pointed out, there are as always multiple ways of achieving this.
Since I wanted to not only output the first principal stress, but actually use it further for my calculations, I wanted to have a DoF vector resembling the field in my mesh. Thus, there are again (at least) 2 options: 
i) introduce a discontinuous field with corresponding DoFs (would require more effort since I am only using continuous FEs)
ii) use a continuous, averaged field, accepting the introduced error (constant per element and averaged over neighbors)

If you want to do ii) as well, it boils down to a classical assembly loop,  and one computes the Cauchy stress tensor (or whatever stress tensor you want to get the principal stresses from) as you did in your momentum balance equation.
This might look something like the following in integration point q:

{
                            Tensor<2,dim> F_s = identity_tensor + gradients_displacement[q];  // displacement gradient
                            Tensor<2,dim> PK2_stress = get_PK2_stress<dim>(F_s); // second Piola-Kirchhoff stress tensor, you define this function based on your constitutive equation
                            Tensor<2,dim> cauchy_stress = symmetrize((1.0/determinant(F_s)) * F_s * PK2_s * transpose(F_s)); // Cauchy stress tensor via Piola transform for finite strain theory

                      const std::array<double,dim> eigenvalues_cauchy = eigenvalues(cauchy_stress);
                            double max_principal_stress = cauchy_stress[0];

}

Optionally, one could always get a global DoF vector via projecting the stress components or the first principal stress (or the quantity you need to plot). This would be much much more accurate than what I suggest in ii).
Going through i) or ii), both ways one ends up with (dis-)continuous global vectors, possibly averaged over neighbors. These are then easily exported via the classical dealii::DataOut<dim> data_out; data_out.add_data_vector(...).
If you want to just output the results, maybe try the dealii::DataPostprocessor<dim> as mentioned by Wolfgang.

Best regards
Richard

[Raúl Aparicio Yuste]

Dear Wolfgang and Richard,

Thank you for your answers. It is good to know how the library works. In my case I am solving an isotropic linear elasticity problem, so I am going to implement it by one of the ways I suggested in my previous comment. Your comments have been really helpful in organising my ideas about Dealii. I will try to come up with the solution and share it to colaborate with this interesting resource. Thank you!

Best

Raúl

[Raúl Aparicio Yuste]

Hello,

Please find attached a script in which I calculate the principal stress components by using the eigen library. This example is in 2D for a particular case, but this piece of code can be modified to generalise the calculation of the principal stress components in 3D. Hope it will be useful. 

Best

Raúl

[Jean-Paul Pelteret]

Hi Raul,

I’m a bit late to the discussion, but maybe this piece of information might still be useful to you. There is a function that computes and returns the eigenvalues and eigenvectors of symmetric tensors: https://dealii.org/developer/doxygen/deal.II/symmetric__tensor_8h.html#aa18a9d623fcd520f022421fd1d6c7a14

 (the eigenvalues of non-symmetric tensors are, in general, complex valued and more complicated to compute, so this has not yet been implemented)

You can do this point-wise and, as previously suggested, could be stored for visualisation using the DataPostprocessTensor class. 

Paraview also has a Tensor glyph (maybe only as an extension?), so maybe that could also be worth investigating.

https://kitware.github.io/paraview-docs/latest/python/paraview.simple.TensorGlyph.html

Best,

Jean-Paul

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/187440c4-5646-43ca-acb3-6ddde9eb46een%40googlegroups.com.
<principal_stress.cpp>

[Raúl Aparicio Yuste]

Hi Jean-Paul,

Thank you for the suggestion, I will take a look in both C++ and Paraview.

Best

Raúl