Hi Jean-Paul,
Thanks for your response. I do think am doing the scaling of the area
correctly, which I based on your answer to my previous question about the
normal vectors. I just did so in a way that is not so clear in the code. To
make it more clear, I have written the code to be
Tensor<1, dim, double> spatial_normal {transpose(invert(deformation_grad))
* material_normal};
spatial_normal /= spatial_normal.norm();
double da_dA {deformation_grad_det * spatial_normal *
transpose(invert(deformation_grad)) * material_normal};
spatial_force[boundary_id - 1] += cauchy * spatial_normal *
fe_face_values.JxW(q_i) * da_dA;
ave_spatial_normal[boundary_id - 1] += spatial_normal;
This gives exactly the same results.
The piola_kirchhoff here is the second Piola-Kirchhoff stress tensor for a
geometrically nonlinear isotropic material. In other words, it's a model
that uses a linear constitutive relation between stress and strain, but
does not make the infinitesimal strain assumption. My understanding of the
second Piola-Kirchhoff stress tensor is that this applies to material
normal vectors to produce pseudo stress/traction vectors that have both a
different orientation and magnitude from the true stress/traction vectors,
as they are relative to the material configuration, and are also in per
unit undeformed area. I have now added in a calculation of the first
Piola-Kichhoff stress tensor (the two-point tensor you mention), the
associated pseudo stress vectors on the boundary faces, and the integrated
normal and shear forces. The magnitude of these forces are equal to the
spatial, or true integrated normal and shear forces, as you said.
const Tensor<2, dim, double> second_piola_kirchhoff {lambda *
trace(green_lagrange_strain) * unit_symmetric_tensor<dim>() + 2 * mu *
green_lagrange_strain};
const Tensor<2, dim, double> first_piola_kirchhoff {deformation_grad *
second_piola_kirchhoff};
second_pseudo_material_force[boundary_id - 1] += second_piola_kirchhoff *
material_normal * fe_face_values.JxW(q_i);
first_pseudo_material_force[boundary_id - 1] += first_piola_kirchhoff *
material_normal * fe_face_values.JxW(q_i);
Still, I do find it unintuitive that the magnitude of the total integrated
force on the boundary is not invariant. As in, why the norm of the summed
true stress vectors is not equal to the norm of the summed pseudo stress
vectors. I thought that by applying the deformation gradient to the first
Piola-Kirchhoff stress tensor you were just reorienting the stress for the
material configuration, but still dealing with the same total forces. I
suppose I will have to reevaluate my intuition about this.
I also still do not know why the analytical result from theory is so
different from the FEM solution with this model, even for such a small
displacement of the boundary (I had a typo in my original post; the beam
theory result gave a shear force with a magnitude of 0.02, while here I
have a shear force magnitude of ~0.15). However, it does not seem the issue
is related to deal.ii, so I will have to consider that further myself.
Thanks so much for your help,
Best,
Alex
On Tuesday, May 11, 2021 at 10:54:31 p.m. UTC+2 Jean-Paul Pelteret wrote:
> Hi Alex,
>
> Yes, at first glance it does look to me like the integration is a problem.
> You are using the same FEFaceValues object to perform the integration in
> both the material and spatial configurations, and this is not quite correct
> with the way that you have things set up. The following (which comes from
> the continuum identity for the volume integral of the divergence of a
> contravariant quantity, i.e. in your case, a stress measure) is equivalent:
>
> F_{t} == \int_{dB_{0}} T_{0} dA = \int_{B_{0}} P . N dA // “Reference
> quantities”, with “P” being the two-point Piola stress tensor
> == \int_{dB_{t}} \sigma . n da = \int_{B_{t}} t da.
> // “Spatial quantities"
>
> That is to say, that the current traction per unit reference area
> integrated over the reference area is equal to the current traction per
> unit current area integrated over the current area. With this operation
>
> spatial_force[boundary_id - 1] += cauchy * spatial_normal *
> fe_face_values.JxW(q_i) *
> deformation_grad_det;
>
> you are integrating the current traction per unit current area over the
> reference area. So you either need to use a mapping and set up
> another FEFaceValues object that would perform integration on the spatial
> configuration, or you can use the following identity for the ratio of
> surface areas to get the correct scaling factor for what is returned by
> fe_face_values.JxW(q_i):
>
> da/dA = det(F) n . F^{-T} . N
>
> That said, there’s also the question about what you’re doing for the
> "material force". Is ” piola_kirchhoff" really the symmetric, fully
> referential Piola-Kirchhoff stress? If so, then no you cannot expect these
> two quantities to match. F_{0} = \int_{dB_{0}} S . N dA is not the same
> force as F_{t}, but is rather a “pseudo-force” that describes some
> reference traction per unit reference area. You need to push forward its
> first index (i.e. compute P = F.S) and then you can get to computing F_{t}
> as stated above.
>
> I hope that helps clear things up a bit.
>
> Best,
> Jean-Paul
>
>
> On 11. May 2021, at 12:46, Alex Cumberworth <alexanderc...@gmail.com>
> wrote:
>
> Hello,
>
> As a test to validate my code, I am solving the equations for
> geometrically nonlinear elasticity (the Saint Venant-Kirchhoff model) for a
> beam with a small displacement boundary condition on the right end such
> that I can compare with Euler-Bernoulli beam theory. I can compare both the
> displacement and the shear force between the FEM solution and the beam
> theory solution. In my FEM integration, I output the normal and shear
> forces for both sides of the beam in both the material and spatial
> reference. The left and right sides are balanced, as expected, but the
> spatial and material forces are not quite equal.
>
> Shouldn't it be the case that spatial and material force is the same? Here
> are the outputted forces for the right side
>
> Right boundary material normal force: 0.0694169
> Right boundary spatial normal force: 0.0724468
>
> Right boundary material shear force: 0.152057
> Right boundary spatial shear force: 0.152864
>
> Further, beam theory gives a shear force with a magnitude of exactly 0.2.
> If I make the displacement smaller the FEM and beam theory shear forces do
> not converge. Is it expected for them to converge?
>
> Below is the deformed system with the stress vectors on the faces
> included. The black grid is the deformed FEM solution, while the solid red
> is the beam theory solution.
>
> <beam.png>
> If there is an issue, I would guess it would be in the integration. In
> converting the material normal vector to the spatial reference, I first
> only applied the inverse transpose of the deformation gradient, and did not
> multiply by the determinant until calculating the force vector. I did this
> so that I can get the unit normal spatial vectors to add up and later
> average so that I have an average normal vector for the whole boundary face
> to calculate the normal and shear force vectors. I have pasted the function
> in below:
>
> template <int dim>
> void SolveRing<dim>::integrate_over_boundaries() {
> QGauss<dim - 1> quadrature_formula(fe.degree + 1);
> FEFaceValues<dim> fe_face_values(
> fe,
> quadrature_formula,
> update_values | update_gradients | update_quadrature_points |
> update_JxW_values | update_normal_vectors);
>
> std::vector<Tensor<1, dim, double>> material_force(2);
> std::vector<Tensor<1, dim, double>> spatial_force(2);
> std::vector<Tensor<1, dim, double>> ave_material_normal(2);
> std::vector<Tensor<1, dim, double>> ave_spatial_normal(2);
> const FEValuesExtractors::Vector displacements(0);
> for (const auto& cell: dof_handler.active_cell_iterators()) {
> for (const auto face_i: GeometryInfo<dim>::face_indices()) {
> const unsigned int boundary_id
> {cell->face(face_i)->boundary_id()};
> if (not(boundary_id == 1 or boundary_id == 2)) {
> continue;
> }
> fe_face_values.reinit(cell, face_i);
> std::vector<Tensor<1, dim, double>> normal_vectors {
> fe_face_values.get_normal_vectors()};
> std::vector<Tensor<2, dim, double>> solution_gradients(
> fe_face_values.n_quadrature_points);
> fe_face_values[displacements].get_function_gradients(
> present_solution, solution_gradients);
> for (const auto q_i:
> fe_face_values.quadrature_point_indices()) {
> const Tensor<2, dim, double> grad_u
> {solution_gradients[q_i]};
> const Tensor<1, dim, double> material_normal {
> normal_vectors[q_i]};
>
> const Tensor<2, dim, double> grad_u_T {transpose(grad_u)};
> const Tensor<2, dim, double> green_lagrange_strain {
> 0.5 * (grad_u + grad_u_T + grad_u_T * grad_u)};
> const Tensor<2, dim, double> piola_kirchhoff {
> lambda * trace(green_lagrange_strain) *
> unit_symmetric_tensor<dim>() +
> 2 * mu * green_lagrange_strain};
>
> ave_material_normal[boundary_id - 1] += material_normal;
> material_force[boundary_id - 1] += piola_kirchhoff *
> material_normal *
> fe_face_values.JxW(q_i);
>
> const Tensor<2, dim, double> deformation_grad {
> grad_u + unit_symmetric_tensor<dim>()};
> const double deformation_grad_det {
> determinant(deformation_grad)};
> const Tensor<2, dim, double> cauchy {
> deformation_grad * piola_kirchhoff *
> transpose(deformation_grad) /
> deformation_grad_det};
>
> Tensor<1, dim, double> spatial_normal {
> transpose(invert(deformation_grad)) *
> material_normal};
> spatial_force[boundary_id - 1] += cauchy * spatial_normal *
> fe_face_values.JxW(q_i) *
> deformation_grad_det;
> spatial_normal /= spatial_normal.norm();
> ave_spatial_normal[boundary_id - 1] += spatial_normal;
> }
> }
> }
> ave_material_normal[0] /= ave_material_normal[0].norm();
> ave_spatial_normal[0] /= ave_spatial_normal[0].norm();
> ave_material_normal[1] /= ave_material_normal[1].norm();
> ave_spatial_normal[1] /= ave_spatial_normal[1].norm();
> const Tensor<1, dim, double> left_material_normal_force = {
> (material_force[0] * ave_material_normal[0]) *
> ave_material_normal[0]};
> const Tensor<1, dim, double> left_material_shear_force = {
> material_force[0] - left_material_normal_force};
> const Tensor<1, dim, double> left_spatial_normal_force = {
> (spatial_force[0] * ave_spatial_normal[0]) *
> ave_spatial_normal[0]};
> const Tensor<1, dim, double> left_spatial_shear_force = {
> spatial_force[0] - left_spatial_normal_force};
> const Tensor<1, dim, double> right_material_normal_force = {
> (material_force[1] * ave_material_normal[1]) *
> ave_material_normal[1]};
> const Tensor<1, dim, double> right_material_shear_force = {
> material_force[1] - right_material_normal_force};
> const Tensor<1, dim, double> right_spatial_normal_force = {
> (spatial_force[1] * ave_spatial_normal[1]) *
> ave_spatial_normal[1]};
> const Tensor<1, dim, double> right_spatial_shear_force = {
> spatial_force[1] - right_spatial_normal_force};
>
> cout << "Left boundary material normal force: "
> << left_material_normal_force.norm() << std::endl;
> cout << "Right boundary material normal force: "
> << right_material_normal_force.norm() << std::endl;
> cout << "Left boundary material shear force: "
> << left_material_shear_force.norm() << std::endl;
> cout << "Right boundary material shear force: "
> << right_material_shear_force.norm() << std::endl;
> cout << "Left boundary spatial normal force: "
> << left_spatial_normal_force.norm() << std::endl;
> cout << "Right boundary spatial normal force: "
> << right_spatial_normal_force.norm() << std::endl;
> cout << "Left boundary spatial shear force: "
> << left_spatial_shear_force.norm() << std::endl;
> cout << "Right boundary spatial shear force: "
> << right_spatial_shear_force.norm() << std::endl;
>
> I will also paste in the core of my assembly loop:
>
> fe_values.reinit(cell);
> cell_matrix = 0;
> cell_rhs = 0;
>
> std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
> cell->get_dof_indices(local_dof_indices);
> const unsigned int n_independent_variables = local_dof_indices.size();
>
> ADHelper ad_helper(n_independent_variables);
> ad_helper.register_dof_values(evaluation_point, local_dof_indices);
> const std::vector<ADNumberType>& dof_values_ad =
> ad_helper.get_sensitive_dof_values();
> ADNumberType energy_ad = ADNumberType(0.0);
>
> std::vector<Tensor<2, dim, ADNumberType>> old_solution_gradients(
> fe_values.n_quadrature_points);
> fe_values[displacements].get_function_gradients_from_local_dof_values(
> dof_values_ad, old_solution_gradients);
>
> for (const unsigned int q_index: fe_values.quadrature_point_indices()) {
> const Tensor<2, dim, ADNumberType> grad_u {
> old_solution_gradients[q_index]};
> const Tensor<2, dim, ADNumberType> grad_u_T {transpose(grad_u)};
> const Tensor<2, dim, ADNumberType> green_lagrange_strain_tensor {
> 0.5 * (grad_u + grad_u_T + grad_u_T * grad_u)};
> ADNumberType t1 = lambda / 2 *
> std::pow(trace(green_lagrange_strain_tensor), 2);
> ADNumberType t2 = mu * double_contract<0, 0, 1, 1>(
> green_lagrange_strain_tensor,
> green_lagrange_strain_tensor);
> ADNumberType pi {t1 + t2};
> energy_ad += pi * fe_values.JxW(q_index);
> }
>
> ad_helper.register_energy_functional(energy_ad);
> present_energy += ad_helper.compute_energy();
> ad_helper.compute_residual(cell_rhs);
> cell_rhs *= -1.0; // RHS = - residual
>
> Any thoughts would be greatly appreciated,
>
> Alex
>
> --
> 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/8b7631a0-5c7e-4bd0-9ff0-79c0d1ed9acdn%40googlegroups.com
>
> <https://groups.google.com/d/msgid/dealii/8b7631a0-5c7e-4bd0-9ff0-79c0d1ed9acdn%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
> <beam.png>
>
>
>
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/f2e4febd-3492-42df-b2b1-ab233f29a8a1n%40googlegroups.com.
