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
Wolfgang Bangerth
On 10/11/25 10:36, Arthur Bawin wrote:
> 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.
want to compare with how ASPECT does it:
https://github.com/geodynamics/aspect/blob/main/source/mesh_deformation/interface.cc
> 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 ?
cases that illustrate what doesn't work, it shouldn't be very difficult
to get those fixed.
> The main, more niche question :
I think you'll have to debug that yourself :-) The usual suggestion
applies as always: make it simple. Try a test case with very low
Reynolds number where flow is stationary, or perhaps even prescribe the
flow. Then you know what the forces will be, and you can think about
what the displacement should be.
Best
W.
Arthur Bawin
I believe I fixed the coupling bug, but I’m a bit perplexed. I narrowed
it down to how the two fields are coupled, and I joined a minimal
example which reproduces the problem I had. It goes as follows :
On the holed domain of the given mesh (cylinder in rectangle), consider
two distinct Poisson equations for two scalar fields u and v, with both
homogeneous Dirichlet BC on the rectangle. The problem is solved in a
monolithic fashion, in the same linear system. On the cylinder, set u
free (zero flux), and say we want the strong coupling v = C*u, without
affecting the values of u.
If the coupling is enforced by simply modifying the rows of the v dofs
on the cylinder, then u is unaffected by the coupling as expected, and
we get v = C*u.
However, if the coupling is enforced with a set of AffineConstraints,
then the v dofs are replaced with C times udofs as the constraints
resolves indirections, which (I believe) modifies the equation for the u
dofs on the cylinder and effectively acts as some kind of penalization.
If the coupling constraint C is large, then u goes to zero on the
cylinder in this example. Is this correct ? And if so, is this the
expected behavior for an AffineConstraints ? I would expect the u dofs
to not be affected by the constraints, acting as « master dofs » (and a
fortiori since is_constrained() for these dofs returns false).
The attached example solves both Poisson problems and compares the
values of u and v on the cylinder (1) when the problem is fully
decoupled (both u and v left free), (2) when the coupling is enforced
with affine constraints, and (3) when enforced by manually adding the
« algebraic » coupling constraints v – C*u = 0.
- In the decoupled case, these values are the same as expected
- In the coupled case with affine constraints, the coupling is satisfied
but the values of u have changed if C != 1 (although we could find other
cases for which they change even for C = 1).
- In the coupled case while manually constraining the system, the
coupling is satisfied and the values of u are unchanged, as required.
I feel like this is straightforward, but I can’t wrap my head around the
fact that naively applying an AffineConstraints for this example
modifies the « source » field.
Thank you for your time,
Arthur
Arthur Bawin
coupling coefficients (./constrained_poisson [1/10/100/...]) illustrates
the discrepancy for the values of u in the "AffineConstraints" case
(that is, the values of u in this case differ more and more from their
expected values in the decoupled case).
Wolfgang Bangerth
>
> If the coupling is enforced by simply modifying the rows of the v dofs
> on the cylinder, then u is unaffected by the coupling as expected, and
> we get v = C*u.
describes, namely
v = C*u,
creates a primary degree of freedom u whose value is what it is, and a
dependent degree of freedom v whose value is then computed from u. But that's
not what the class does. It represents constraints that you should think have
the form
v - C*u = 0
in which all degrees of freedom are logically considered having the same category.
I put this into words here:
https://github.com/dealii/dealii/pull/18945
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/