Deal deal.ii community,

I am writing a monolithic fluid-structure interaction test case in an 
arbitrary Lagrangian-Eulerian (ALE) formulation with the toolbox, and I 
have a few interrogations.

The setting is the following : I consider an incompressible flow past a 
massless and rigid cylinder, which is attached to a spring of stiffness *k*. 
The cylinder dynamics is then simply a balance between the spring force F_s 
and the force exerted by the fluid F_f. The fluid exerts a force on the 
cylinder, which moves and moves the mesh with it, which in turn modifies 
the flow geometry, and so on. To work in a fully implicit setting, the 
no-slip condition on the cylinder is enforced with a Lagrange multiplier, 
allowing to recover the fluid forces on the cylinder which are up to sign 
the integral of the Lagrange multiplier. To move the mesh, I use a 
pseudo-solid analogy, where the mesh movement is obtained by solving the 
linear elasticity equations with arbitrary and equal Lamé coefficients (say 
1 for both). The position boundary condition on the cylinder is given by 
the coupling between position and Lagrange multiplier DOFs, which is the 
equation F_s + F_f = 0.

The problem unknowns are thus the fluid velocity u, the pressure p, the 
mesh position x and the Lagrange multiplier lambda, the latter defined only 
on the boundary of the cylinder. The elasticity equation is solved in a 
Lagrangian setting on the initial mesh, using a MappingFE with linear 
FE_SimplexP, whereas the incompressible Navier-Stokes equations are solved 
on the moving mesh, using a MappingFEField attached to the position *x* 
part of the solution vector.

I am using triangles/tetrahedra, Taylor-Hood P2-P1 elements for the 
velocity-pressure pair, a linear (P1) representation of the geometry x and 
quadratic (P2) Lagrange multiplier. A similar setting (although with P2 
geometry) is described here : https://doi.org/10.1016/j.jcp.2020.109734 , 
with the difference that mesh movement is treated by interpolating the mesh 
velocity instead of using the pseudo-solid approach.

A first general question :

- The Lagrange multiplier should be defined only on the cylinder boundary, 
a codimension 1 space. What is the simplest way of dealing with a <dim-1> 
finite element space within a FESystem that contains the other spaces of 
dimension dim ? Currently, and to use a single FESystem and DOFHandler, I 
define lambda as a P2 field over the whole domain and constrain all of its 
non-cylinder DOFs to zero, which works but is not very elegant. I believe 
that FENothing could be used to this end, but when I tried it seems that it 
is not fully implemented for simplices, is that correct ?

The main, more niche question :

- The fluid-structure coupling on the cylinder is as follows : the lambda 
DOFs are computed to enforce *u - w = 0* weakly on the cylinder where *w* 
is the mesh velocity, *w* is computed from the position DOFs using the same 
BDF formula as for *dudt*, and the position DOFs are constrained to the 
lambda DOFs to represent the cylinder dynamics, using an AffineConstraint. 
This constraint is x_c = X_c – 1/k * int_c lambda dl, where X_c is the 
initial position on the cylinder. The problem I currently have is that the 
weak no-slip is not enforced as soon as I apply the coupling between x and 
lambda : the affine constraint is satisfied (I check that the prescribed 
displacement is indeed -1/k times the integral of lambda), but *u != w*. If 
the position is not coupled, for instance if the movement of the cylinder 
is zero or prescribed (e.g., a sine function), then the no-slip is enforced 
exactly. What is very weird is that when the coupling is enabled, the (L2 
and Linf) error on the no-slip scales with the magnitude of the Lamé 
parameters of the mesh pseudo-solid, although these should not interact.

My initial Newton guess for the position is x = X, which satisfies the 
inhomogeneous part of the constraint, so I only enforce the homogeneous 
part in the Newton resolution by using 

zero_constraints.distribute_local_to_global(local_matrix,

local_dof_indices,

system_matrix);

for the subsequent iterations and time steps, and similarly for the rhs, 
not using the call with 5/6 arguments.

Here are the functions which compute the coupling coefficients and create 
the position-lambda constraints :

template <int dim>
  void FSI<dim>::create_position_lambda_coupling_coefficients(
    const unsigned int boundary_id)
  {
    // l_lower is the lower lambda component in the FESystem which contains
    // u, p, x, lambda, i.e., 2 * dim + 1.
    const FEValuesExtractors::Vector lambda(l_lower);

    //
    // Compute the weights c_ij.
    // Done only once as cylinder is rigid and those weights will not 
change.
    //
    std::vector<std::map<types::global_dof_index, double>> coeffs(dim);

    FEFaceValues<dim> fe_face_values_fixed(*fixed_mapping,
                                           fe,
                                           face_quadrature,
                                           update_values | 
update_JxW_values);

    std::vector<types::global_dof_index> face_dofs(fe.n_dofs_per_face());

    for (const auto &cell : dof_handler.active_cell_iterators())
    {
      if (!cell->is_locally_owned())
        continue;
      for (const auto i_face : cell->face_indices())
      {
        const auto &face = cell->face(i_face);
        if (!(face->at_boundary() && face->boundary_id() == boundary_id))
          continue;

        fe_face_values_fixed.reinit(cell, face);
        face->get_dof_indices(face_dofs);

        for (unsigned int q = 0; q < face_quadrature.size(); ++q)
        {
          const double JxW = fe_face_values_fixed.JxW(q);

          for (unsigned int i_dof = 0; i_dof < fe.n_dofs_per_face(); 
++i_dof)
          {
            const unsigned int comp =
              fe.face_system_to_component_index(i_dof, i_face).first;

            // Account for ghost DoF (not only owned), which contribute to 
the
            // integral on this element
            if (!locally_relevant_dofs.is_element(face_dofs[i_dof]))
              continue;

            if (is_lambda(comp))
            {
              const types::global_dof_index lambda_dof = face_dofs[i_dof];

              // Get cell dof index from this face dof index
              const unsigned int i_cell_dof = fe.face_to_cell_index(i_dof, 
i_face);
              const unsigned int d     = comp - l_lower;
              const double       phi_i = 
fe_face_values_fixed.shape_value(i_cell_dof, q);

              if(coeffs[d].count(lambda_dof) == 0)
                coeffs[d][lambda_dof] = 0.;

              coeffs[d][lambda_dof] += - phi_i * JxW / 
param.spring_constant;
            }
          }
        }
      }
    }

    // Move map coefficients to vector member
    position_lambda_coeffs.resize(dim);
    for (unsigned int d = 0; d < dim; ++d)
    {
      position_lambda_coeffs[d].insert(position_lambda_coeffs[d].end(),
                                       coeffs[d].begin(),
                                       coeffs[d].end());
    }
  }

  template <int dim>
  void FSI<dim>::create_position_lambda_constraints(const unsigned int 
boundary_id)
  {
    position_constraints.clear();
    position_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);

    std::vector<types::global_dof_index> face_dofs(fe.n_dofs_per_face());

    for (const auto &cell : dof_handler.active_cell_iterators())
    {
      if (!cell->is_locally_owned())
        continue;
      for (const auto i_face : cell->face_indices())
      {
        const auto &face = cell->face(i_face);
        if (!(face->at_boundary() && face->boundary_id() == boundary_id))
          continue;

        face->get_dof_indices(face_dofs);

        for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
        {
          if (!locally_owned_dofs.is_element(face_dofs[i]))
            continue;

          const unsigned int comp =
            fe.face_system_to_component_index(i, i_face).first;

          if (is_position(comp))
          {
            const unsigned int d = comp - x_lower;
            position_constraints.add_line(face_dofs[i]);
            position_constraints.add_entries(face_dofs[i], 
position_lambda_coeffs[d]);
          }
        }
      }
    }
    position_constraints.close();
  } 

Sorry for the very long post, and thank you all for your time.

Arthur

-- 
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 [email protected].
To view this discussion visit 
https://groups.google.com/d/msgid/dealii/0e7fdead-836e-4cfc-b025-7792f095af38n%40googlegroups.com.

Reply via email to