Pressure load on a plate upper surface in solid mechanics

[Léonhard YU]

Dear all,

 I am working on impose constant pressure load on a plate's upper surface. The code is derived from step-18 and attached below. But I found that the diplacement is not right, and it varies when the number of face elements changes. I am thinking about something wrong with the loading but I do not have any idea. Thanks a lot if you can help me.

The constant pressure is 1e-2 and the boundary_id for loading is bc_uppersurface. The codes attached below:

for (const auto &cell : dof_handler.active_cell_iterators())
    if (cell->is_locally_owned())
    {
        cell_matrix = 0;
        cell_rhs    = 0;

        fe_values.reinit(cell);
        for (unsigned int i = 0; i < dofs_per_cell; ++i)
        for (unsigned int j = 0; j < dofs_per_cell; ++j)
            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
                const SymmetricTensor<2, dim>
                eps_phi_i = get_strain(fe_values, i, q_point),
                eps_phi_j = get_strain(fe_values, j, q_point);

                cell_matrix(i, j) += (eps_phi_i *            //
                                    stress_strain_tensor * //
                                    eps_phi_j              //
                                    ) *                    //
                                    fe_values.JxW(q_point); //
            }

    // pressure load
        for (const auto &face : cell->face_iterators())
        {
            if (face->at_boundary() && (face->boundary_id() == bc_uppersurface))
            {
            fe_face_values.reinit(cell, face);
            for (unsigned int i = 0; i < dofs_per_cell; ++i)
            {
                const unsigned int component_i = fe.system_to_component_index(i).first;
                for (unsigned int face_q_point = 0; face_q_point < face_n_q_point; ++face_q_point)
                {
                Tensor<1,dim> dir;
                dir = fe_face_values.normal_vector(face_q_point);
                Point<dim> point_dof = fe_face_values.quadrature_point(face_q_point);
                const double magnitude = 1e-2;
                const Tensor<1, dim> traction = -1.0 * magnitude * dir;
                cell_rhs(i) +=
                            fe_face_values.shape_value(i, face_q_point) * traction[component_i] *
                        fe_face_values.JxW(face_q_point);
                }
            }
            }          
        }
    }

Best Regards

[陈敏]

Hi  Léonhard

Your solution may not Converge to a mesh-independent solution.

Best 

Chen

Léonhard YU <leonh...@gmail.com> 于2024年6月8日周六 18:14写道：

--
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/d9ed8d3d-aa19-4367-9505-ca9d9656d45cn%40googlegroups.com.

[陈敏]

One more thing: 

you should loop over q point before loop over dof.

Léonhard YU <leonh...@gmail.com> 于2024年6月8日周六 18:14写道：

Dear all,

--

[Wolfgang Bangerth]

On 6/8/24 04:14, Léonhard YU wrote:
>  I am working on impose constant pressure load on a plate's upper
> surface. The code is derived from step-18 and attached below. But I
> found that the diplacement is not right, and it varies when the number
> of face elements changes. I am thinking about something wrong with the
> loading but I do not have any idea. Thanks a lot if you can help me.

Leonhart:
the code looks not unreasonable to me (though I don't know what
get_strain() does). But as a general rule, when you're trying to debug
programs, it's very useful to be specific about what "the diplacement is

not right, and it varies when the number of face elements changes"

actually means. That's true when describing the issue to others, but
also to yourself. So:
* How exactly is it wrong?
* What is it you are expecting?
* What is it you are seeing?
* What is the difference?
* Do you know things that the solution should satisfy, such as that it
needs to be left-right symmetric because of the chosen forcing terms and
geometry?
* Does the computed solution satisfy these condition?

These are all questions that will lead you to the place where your
program is doing something wrong.

Best
W.