Integrated material and spatial traction forces on boundary not equal

[Alex Cumberworth]

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.

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

[Jean-Paul Pelteret]

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.

--
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.
<beam.png>

[Alex Cumberworth]

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

[Michael Lee]

Hi Alex,

I'm learning deal.ii and trying do the similar verification. If it is possible for you to share the code with me?

Thank you!

Michael

[Alex Cumberworth]

Hello,

I have attached the most recent version of my code here. I have tried to make setting boundary conditions in the parameter file more convenient for myself; you can set boundary domains, boundary conditions that use these domains, and stages if the displacement is large. However, there are not many comments, so you may want to just remove this part for your own purposes.

However, I have been a bit surprised in my comparisons at how fine a mesh is required to achieve convergence with the beam theory result. I am now using a beam that is 1 x 2 x 20, and using the subdivided rectangle helper function, I set the number of subdivisions to be 5 and 10 for the width and height, respectively. I then varied the number of subdivision in the length between 10 and 1000. The beam theory result is that the shear force has a magnitude of 0.001 for a displacement on the right side of 0.1. Even at 1000 subdivisions, the FEM result is 0.00113 (from 0.00129 at 500). The system has 200 000 degrees of freedom, and the result is still off by 13%. Is it expected that even to calculate a shear force in this simple problem that such a large number of degrees of freedom are required?

Best,

Alex

[Michael Li]

Alex,

 

Thank you for sharing your codes. I have some compiling errors relating to “EnergyFunctional’:

solve_ring_nonlinear.cpp:520:43: error: ‘EnergyFunctional’ in namespace ‘dealii::Differentiation::AD’ does not name a template type

  520 |     using ADHelper = Differentiation::AD::EnergyFunctional<

      |                                           ^~~~~~~~~~~~~~~~

Can you help me with that?

 

For the large error, I noticed you’re using linear element. I encountered the same large error when comparing it with Abaqus with FE_Q(1). But the error came down with less grids when I used higher finite element FE_Q<dim>(2). I remember the deflection of beam is a cubic function of coordinate. You may try that see if it improves.

 

Best,

Michael

 

From: Alex Cumberworth
Sent: Wednesday, June 9, 2021 9:54 AM
To: deal.II User Group
Subject: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

 

Hello,

 

I have attached the most recent version of my code here. I have tried to make setting boundary conditions in the parameter file more convenient for myself; you can set boundary domains, boundary conditions that use these domains, and stages if the displacement is large. However, there are not many comments, so you may want to just remove this part for your own purposes.

 

However, I have been a bit surprised in my comparisons at how fine a mesh is required to achieve convergence with the beam theory result. I am now using a beam that is 1 x 2 x 20, and using the subdivided rectangle helper function, I set the number of subdivisions to be 5 and 10 for the width and height, respectively. I then varied the number of subdivision in the length between 10 and 1000. The beam theory result is that the shear force has a magnitude of 0.001 for a displacement on the right side of 0.1. Even at 1000 subdivisions, the FEM result is 0.00113 (from 0.00129 at 500). The system has 200 000 degrees of freedom, and the result is still off by 13%. Is it expected that even to calculate a shear force in this simple problem that such a large number of degrees of freedom are required?

 

 

Best,

Alex

 

On Monday, June 7, 2021 at 3:39:44 p.m. UTC+2 lian...@gmail.com wrote:

Hi Alex,

I'm learning deal.ii and trying do the similar verification. If it is possible for you to share the code with me?

 

Thank you!

Michael

On Tuesday, May 11, 2021 at 4:46:55 AM UTC-6 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.

--

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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/K-lMxbtZUdQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4cade24b-48aa-432b-a590-a75535c0f2f2n%40googlegroups.com.

 

[Alex Cumberworth]

Hi Michael,

I think that is the first point in my code that something from the automatic differentiation module is referred to, which is a part of trilinos. Was your version of deal.ii compiled with trilinos, and if so, was your version of trilinos compiled with sacado?

To answer the question I posed, I did some more tests and found that refining in the width and height of the beam has no effect on the convergence of the shear force. I went down to just a single division in z and 2 in y (in other words there are just two faces when looking down the end of the beam), and could get the same level of agreement with beam theory with one order of magnitude fewer degrees of freedom. I guess I am still mildly surprised how refined I had to made the x direction, but perhaps if I also used adaptive refinement things would further improve.

Best,

Alex

[Wolfgang Bangerth]

On 6/11/21 4:04 AM, Alex Cumberworth wrote:
>
> To answer the question I posed, I did some more tests and found that refining
> in the width and height of the beam has no effect on the convergence of the
> shear force. I went down to just a single division in z and 2 in y (in other
> words there are just two faces when looking down the end of the beam), and
> could get the same level of agreement with beam theory with one order of
> magnitude fewer degrees of freedom. I guess I am still mildly surprised how
> refined I had to made the x direction, but perhaps if I also used adaptive
> refinement things would further improve.

You might also want to try quadratic or cubic elements. For smooth solutions,
that often does wonders!

Best
W>


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

[Alex Cumberworth]

Hi Wolfgang,

I tried a quadratic element, and was able to reduce the number of degrees of freedom for the same level of convergence by another order of magnitude. Thanks!

Best,

Alex

[Michael Lee]

Hi Alex,

I just got a chance to have Trilinos installed but the same error occurs regarding EnergyFunctional. I got the following notice when trying to test step-31:

Error! This tutorial requires a deal.II library that was configured with

  the following options:

      DEAL_II_WITH_TRILINOS = ON

  However, the deal.II library found at /usr/local was configured with these

  options

      DEAL_II_WITH_TRILINOS = OFF

  which conflict with the requirements.

-- Configuring incomplete, errors occurred!  

Is there a way to let deal.II know Trillinos is installed? Do I have to reinstall deal.ii?

Or, can I compile the code without using Trillinos with some modification?

Thanks!

Michael

[Wolfgang Bangerth]

On 6/16/21 4:07 PM, Michael Lee wrote:
>
> Is there a way to let deal.II know Trillinos is installed? Do I have to
> reinstall deal.ii?
> Or, can I compile the code without using Trillinos with some modification?

No, you have to re-install deal.II so that at compile time, it knows that
Trilinos exists.

Best
W.

[Michael Li]

Thanks, W. for confirming that! I reinstalled deal.II after installing Trillinos. Step 31 runs well.

--

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 a topic in the Google Groups "deal.II User Group" group.

To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/K-lMxbtZUdQ/unsubscribe.

To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/06a30833-0914-ed13-b4e1-dfd9c2c3b72d%40colostate.edu.