How to enforce fluid and solid dynamic coupling in fluid-structure interactions using the finite element method?

Question

I apologize in advance if the question has been posted before or if it sounds a bit naive.

I am writing my own code in MATLAB for a staggered finite element solver for fluid-structure interaction problems with small displacements. So far, I have benchmarked the fluid and solid solvers separately and both are working fine.

At this point I am stuck in how to couple the forces between the fluid and the solid domains. I know that the correct approach is to balance the stresses on the interface from both the solid and the fluid domains as follows:

\begin{equation} \boldsymbol{\sigma}_s \cdot \mathbf{n}_s = \boldsymbol{\sigma}_f \cdot \mathbf{n}_f \end{equation}

But so far I am unable to translate that condition into the discretized weak form and later into code.

To benchmark my force coupling, I am using two set of values from COMSOL Multiphysics (same problem setup): i. Compare the total reaction on the solid (in $x$ and $y$ directions). ii. Compare the maximum displacement after solving the solid system.

To find an easier solution for the coupling, I tried the following two approaches:

1. After solving the fluid system for the fluid state variables, I used the fluid finite element system of equations directly to calculate the forces (see below, Reddy and Gartling 2010 P166). I picked the negative of those forces that correspond to the solid interface nodes.

I understand that the fluid uses Eulerian description while the solid uses the Lagrangian. But since I am assuming infinitesimal deformations, the stress tensors should not be that different.

This approach gives less than 2% error in total reactions but between 5%-10% in maximum displacement.

2. Using the traction boundary condition formulation on the fluid to calculate the forces on the fluid-structure interface (see below). I picked the negative of these values that correspond to the interface nodes, but I am still getting wrong total reactions and maximum displacement.

In formulating this condition, I followed Reddy and Gartling 2010 P184.

My questions are:

a. Are approaches 1 and 2 theoretically correct to achieve the dynamic coupling between the fluid and solid domains?

b. Are there any resources on how to translate the stress balance condition at the interface into the discretized weak form and/or into code?

Answer 1

You could alternatively start off with a matching discretization at the fluid-solid interface and a monolithic scheme. While writing down the variational formulation of your coupled fluid-solid problem, the fluid and stress tractions can be then made to automatically cancel each other. This is possible by an appropriate choice of the test functions for the fluid and solid momentum balance equations. see the Chapter 5 of this detailed book or paper for more details.

Answer 2

Since I managed to solve the problem I had with force coupling in fluid-structure interactions, I thought it would be best to post a detailed answer for anyone facing the same problem.

As in any textbooks on the subject, the equation to couple the forces between the fluid and structure in the strong form is:

\begin{equation} \boldsymbol{\sigma}_s \cdot \mathbf{n}_s = \boldsymbol{\sigma}_f \cdot \mathbf{n}_f \end{equation}

This equation must be transformed into the weak form, then discretized in order be applied in a finite element solver.

The two approaches I was considering to solve this problem were:

1. After solving the fluid system for the fluid state variables, I used the fluid finite element system of equations directly to calculate the forces (see below, Reddy and Gartling 2010 P166). I picked the negative of those forces that correspond to the solid interface nodes. This method doesn't have any theoretical foundation to back it (more of a hunch), so I am not surprised it doesn't work, and I didn't spend too much time figuring out why it doesn't work.

2. Using the traction boundary condition formulation on the fluid to calculate the forces on the fluid-structure interface (see below). I picked the negative of these values that correspond to the interface nodes, but I am still getting wrong total reactions and maximum displacement. This approach is in fact the correct one, the errors I was getting were unfortunately due to a mistake in the textbook I was following (Reddy and Gartling 2010). This particular mistake in the textbook is in Eq. 4.5.22b in page 184 in the 2010 edition, and hence it propagated to the above two expressions for the $x$ and $y$ forces (Eqs. 4.5.23a and 4.5.23b in the textbook). The issue is that the authors used shape functions to calculate the normals on the finite element edge of interest. However, using shape functions doesn't give consistent directions for the normals. In other words, the normal calculated this way on one particular edge would point inwards, but on a different edge of the same finite element would point outwards. To overcome this issue, I had to go back to simple geometry concepts to calculate normals from nodal coordinates while making sure consistent nodal ordering is used for each element. See the below hand-written note on how to I did it.

To answer the questions I posted at the end of my original post:

a. Are approaches 1 and 2 theoretically correct to achieve the dynamic coupling between the fluid and solid domains?

Approach one doesn't have a theoretical foundation, so it's most probably theoretically incorrect. Approach 2 is the theoretically and practically correct one, I just had the issue with the wrong equations to calculate the normals as detailed above.

b. Are there any resources on how to translate the stress balance condition at the interface into the discretized weak form and/or into code?

In Reddy and Gartling 2010 page 184, I used Eqs. 4.5.21a and 4.5.21b to calculate the forces in the $x$ and $y$ directions respectively. Again, don't use Eq. 4.5.22b to calculate the normals as it's incorrect.

There is one more issue that is relevant to this topic; that is the discontinuity of stresses across element edges in the classical Galerkin finite element method with Lagrangian shape functions. Since the fluidic stresses depend on the first derivatives of the velocity components and continuity is only imposed on the primitive variables (velocity and pressures) themselves not their first derivatives, there would be an inaccuracy in the calculation of the stresses (or fluxes in other physics). There is two approaches to overcome this issue; either impose the coupling condition in the weak form, or use a smoothing technique to approximate the stresses/fluxes across finite element faces/edges. I haven't yet implemented either of those methods, but I am working on them at the moment. There is also some information in COMSOL Multiphysics Reference Manual under "Calculating Accurate Fluxes" if you are interested in reading more on this issue.

Finally, I wrote a manuscript on force coupling in topology optimization of fluid-structure interactions but unfortunately it was rejected because of this stress discontuity issue. It's available as a preprint on arxiv here if anyone is interested. At least, the force coupling derivation is correct in this manuscript but like I said before you have to either live with the stress discontinuity or figure out a way to implement/approximate accurate fluxes.