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
El lunes, 21 de noviembre de 2022 a las 8:41:38 UTC+1, Raúl Aparicio Yuste
escribió:
> 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
>
> On Saturday, November 19, 2022 at 5:43:14 PM UTC+1 [email protected]
> wrote:
>
>> 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
>>
>>
>> [email protected] schrieb am Donnerstag, 17. November 2022 um 16:21:30
>> UTC+1:
>>
>>> 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
>>>
>>
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/187440c4-5646-43ca-acb3-6ddde9eb46een%40googlegroups.com.
#include <eigen3/Eigen/Dense> // library to compute eigenvalues
using Eigen::MatrixXd;
using Eigen::EigenSolver;
using Eigen::VectorXcd;
MatrixXd A(2,2); // matrix declaration
// Stress tensor in 2D
A(0,0) = sigma_x;
A(0,1) = tau_xy;
A(1,0) = tau_xy;
A(1,1) = sigma_y;
// Compute eigen values
EigenSolver<MatrixXd> es(A);
auto eigenvalues_A = es.eigenvalues();
auto eigenvectors_A = es.eigenvectors();
// Consider the first eigenvalue:
double lambda_1A = es.eigenvalues()[0].real();
// Consider the second eigenvalue:
double lambda_2A = es.eigenvalues()[1].real();
// Consider the first eigenvector:
VectorXcd v_1 = es.eigenvectors().col(0);
// Consider the second eigenvector:
VectorXcd v_2 = es.eigenvectors().col(1);
// Initialize variables
double max_principal_stress = 0.0;
double min_principal_stress = 0.0;
double max_principal_stress_direction_1 = 0.0;
double max_principal_stress_direction_2 = 0.0;
double min_principal_stress_direction_1 = 0.0;
double min_principal_stress_direction_2 = 0.0;
// organize eigen values and eigen vectors, so sigma_I > sigma_II
if (lambda_1A >= lambda_2A) {
max_principal_stress = lambda_1A;
min_principal_stress = lambda_2A;
max_principal_stress_direction_1 = v_1(0).real();
max_principal_stress_direction_2 = v_1(1).real();
min_principal_stress_direction_1 = v_2(0).real();
min_principal_stress_direction_2 = v_2(1).real();
} else if (lambda_1A < lambda_2A) {
max_principal_stress = lambda_2A;
min_principal_stress = lambda_1A;
max_principal_stress_direction_1 = v_2(0).real();
max_principal_stress_direction_2 = v_2(1).real();
min_principal_stress_direction_1 = v_1(0).real();
min_principal_stress_direction_2 = v_1(1).real();
}
// Project principal components in order to plot principal stresses in X,Y coordinates
// Make use of scalar product and vector projection in 2D
double stress_I_projection_x = 0.0;
double stress_I_projection_y = 0.0;
double stress_II_projection_x = 0.0;
double stress_II_projection_y = 0.0;
stress_I_projection_x = max_principal_stress * max_principal_stress_direction_1; // stress_I vector projected in X axis
stress_I_projection_y = max_principal_stress * max_principal_stress_direction_2; // stress_I vector projected in Y axis
stress_II_projection_x = min_principal_stress * min_principal_stress_direction_1; // stress_II vector projected in X axis
stress_II_projection_y = min_principal_stress * min_principal_stress_direction_2; // stress_II vector projected in Y axis
