Coupling issue in fluid-structure interaction test
Arthur Bawin
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