Boundary condition from nodal traction values

[Lucky Qin]

Hello everyone,

I would like to apply a traction boundary condition on the surface of a cantilever beam. Suppose I have a vector named traction[dof_handler.n_dofs()], which gives the various traction on each dof. I implemented the boundary condition as follows:

            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 = 0; face < GeometryInfo<dim>::faces_per_cell; ++face )
                {
                    if ( cell->face( face )->at_boundary() == true && cell->face( face )->boundary_id() != 0 )   //boundary_id() == 0  means the Dirichlet boundary condition
                    {
                        fe_face_values.reinit( cell, face );

                        std::vector<unsigned int> local_face_dof_indices( dofs_per_face );
                        cell->face( face )->get_dof_indices( local_face_dof_indices );


                        assert( dofs_per_face == fe_face_values.n_quadrature_points * dim );


                        for ( unsigned int q = 0; q < fe_face_values.n_quadrature_points; ++q )
                        {
                            cell_rhs( i ) += traction( local_face_dof_indices[q * dim + component_i] ) * fe_face_values.shape_value( i, q ) * fe_face_values.JxW( q );
                        }
                    }
                }
            }

At present, I set the traction to be nonzero in only one of the surfaces of a cantilever beam to make quantitative comparison with the previous literature. However, the results do not seem to be correct.

Another thing I would like to know is what this traction may mean in this code. Is this pressure or force?

I would like to ask for some help from you about these two questions.

Best regards，

Qin

[Jean-Paul Pelteret]

Dear Qin,

If I understand this statement 

 I have a vector named traction[dof_handler.n_dofs()], which gives the various traction on each dof

correctly, you have a pre-assembled global traction vector (of size n_dofs) that you wish to add to the system. Why is it that you want to “assemble” this vector again? Can you not just call system_rhs += traction_vector? If the Neumann boundaries overlapped with the Dirichlet boundaries, then you might be able to first call constraints.set_zero(traction_vector) so zero out the contributions from the overlapping DoFs.

If that’s not a viable solution then one problem in the code that you’ve shown is that your assembly itself is incorrect. You need to extract the traction vector at the face quadrature point and assemble with that. Currently you’re performing assembly at each face quadrature point using the value of some element of the total traction vector as located at a support point. Remember that if you were to integrate some local traction vectors over the Neumann boundaries then you should end up with the global traction vector that you have precomputed. You can get the quadrature point values of the local traction vector from the global vector using FEValues::get_function_values() or a similar call.

I hope that this helps.

Best,

Jean-Paul

--
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.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 10/07/2018 09:29 AM, Lucky Qin wrote:
>
> I would like to apply a traction boundary condition on the surface of a
> cantilever beam. Suppose I have a vector named

> /traction[dof_handler.n_dofs()/*], *which gives the various traction on
> each dof.

To augment the answer Jean-Paul has already given:

The way you see things when you define such a traction vector is not the
way you should be thinking of tractions in the finite element context.
There, you should think of a traction boundary condition as a function
g(x) so that

n . sigma(u(x)) = g(x)

where sigma(u) gives you the stress at a point u. In linear elasticity,
for example, this would be something like

sigma(u) = lambda eps(u) + 2*mu*div(u)*I

if I recall correctly. The point I want to make is that the traction is
defined *for every point of the boundary*. This then allows you to
compute contributions to the right hand side of the discrete problem by
approximating the boundary integral

T_i = \int_{\partial \Omega} g(x) phi_i(x) dx

via quadrature. For this to work, you have to have the traction at every
point of the boundary.

On the other hand, you only provide the traction at individual points
(assuming that's how you have defined your 'traction' array), which is
not good enough to actually do the integration.

The reason for all of this is that you will want to work on a sequence
of successively finer meshes to explore convergence of your solution.
For this, it is not enough to just know the traction at individual
points -- you need to be able to evaluate it everywhere.

Best
W.

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

[Lucky Qin]

Dear Jean-Paul,

Thank you for the answer. Sorry for the late response.

I tried to use 'system_rhs += traction_vector' and it works! The code is much easier now.

I do have a pre-assembled global traction vector. This vector is obtained from a fluid solver (based on Lattice Boltzmann method). Since the fluid solver is a discrete solver. I only have traction of the boundaries on vertices.

Thank you again for the suggestions!

Best,

Qin

在 2018年10月8日星期一 UTC+8下午3:51:24，Jean-Paul Pelteret写道：

[Lucky Qin]

Dear Dr. Bangerth,

Thank you for the explanation. I am trying to study fluid-structure interaction and use Dealii as  the solid solver for the problem. However, I only have a traction vector at individual points. This is because the traction vector is obtained from a discrete fluid solver (Lattice Boltzmann solver). The traction vector is discontinuous and I do not have a function to describe it.

Thanks,

Qin

在 2018年10月9日星期二 UTC+8上午3:46:36，Wolfgang Bangerth写道：