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 <alexandercumberwo...@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/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <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
> <mailto:dealii+unsubscr...@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/9685B0E3-E98E-4717-93B2-1289BDA83364%40gmail.com.
