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.
[image: 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+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.
