Finite Element Discretization for Fluid-Structure Interaction
A flexible valve leaf bends under blood flow. A bridge cable vibrates in wind. A hydrofoil twists as water loads rise. In each case, the fluid changes the structure, and the structure changes the fluid. That two-way exchange is the core difficulty of fluid-structure interaction, or FSI.
The finite element method is a natural tool for this problem because it can represent complex shapes, moving boundaries, and different material laws in one framework. Yet FSI is not just “fluid FEM plus solid FEM.” The interface couples two physics models that use different descriptions of motion, different stability needs, and often different time scales.
This post walks through Finite Element Discretization for Fluid-Structure Interaction from the modeling level down to the algebraic system. The focus is on the practical choices that shape a working method: domains, weak forms, interface conditions, element spaces, time integration, and coupling strategy.
FSI starts with two descriptions of motion
The first challenge in FSI is that the fluid and structure are usually described in different ways.
A solid is often modeled in a Lagrangian frame. Material points are tracked as they move. If a point starts at a reference coordinate \(X\), its current position is
\[
x = X + d_s(X,t)
\]
where \(d_s\) is the structural displacement.
A fluid is often modeled in an Eulerian frame. The mesh sits in space while fluid moves through it. Velocity and pressure are fields such as
\[
u_f(x,t), \quad p_f(x,t)
\]
That setup works well for fixed domains. FSI breaks the comfort zone because the fluid boundary moves with the structure.
A common finite element treatment uses an Arbitrary Lagrangian-Eulerian formulation, usually called ALE. The fluid mesh moves, but not necessarily with the fluid particles. It follows the interface exactly and deforms inside the fluid domain to preserve element quality.
In ALE form, the convective term uses the relative velocity between the fluid and the moving mesh:
\[
(u_f - w)\cdot \nabla u_f
\]
where \(w\) is the mesh velocity.
This small change matters. If the mesh follows the boundary but the equations ignore mesh motion, the discretization no longer represents the intended physics.
The interface carries the main coupling
Let \(\Gamma_{fs}\) be the shared boundary between the fluid and the solid. Two key conditions hold there.
The kinematic condition says the fluid velocity matches the structural velocity:
\[
u_f = \frac{\partial d_s}{\partial t}
\quad \text{on } \Gamma_{fs}
\]
For viscous flow, this is the no-slip and no-penetration condition on a moving wall.
The dynamic condition says tractions balance:
\[
\sigma_f n_f + \sigma_s n_s = 0
\quad \text{on } \Gamma_{fs}
\]
Here, \(\sigma_f\) is the fluid stress, \(\sigma_s\) is the solid stress, and the normals point outward from their respective domains.
These two conditions create the feedback loop. The structure moves the fluid boundary. The fluid stress loads the structure. A finite element discretization must preserve that exchange as accurately and stably as possible.
The weak form sets the stage for finite elements
Finite element methods start by writing the governing equations in weak form. Instead of demanding that the differential equation hold point by point, the weak form asks it to hold against a set of test functions. This allows lower regularity, handles complex geometries, and makes boundary terms explicit.
For an incompressible Newtonian fluid in a moving ALE domain, a common model is
\[
\rho_f \left(
\frac{\partial u_f}{\partial t}\bigg|_{\chi}
+ ((u_f-w)\cdot\nabla)u_f
\right)
\nabla\cdot\sigma_f = f_f
\]
\[
\nabla\cdot u_f = 0
\]
with stress
\[
\sigma_f = -p_f I + 2\mu_f \varepsilon(u_f)
\]
and strain-rate tensor
\[
\varepsilon(u_f)=\frac{1}{2}(\nabla u_f+\nabla u_f^T)
\]
For the solid, the exact equation depends on the material model and whether deformations are small or large. A simple linear elastic model in reference coordinates is
\[
\rho_s \frac{\partial^2 d_s}{\partial t^2}
\nabla\cdot\sigma_s(d_s) = f_s
\]
For large deformation, the solid is usually written in a total Lagrangian or updated Lagrangian form with deformation gradient \(F\), Green-Lagrange strain, and a chosen constitutive law.
The weak form appears after multiplying by test functions and integrating over each domain. In the fluid momentum equation, integration by parts moves derivatives from stress to the test function and exposes boundary traction terms. On the interface, those traction terms meet the structural traction terms.
That is where FSI differs from two separate simulations. The boundary term is not just a load or a constraint. It is part of a coupled variational statement.
Matching and nonmatching meshes
At the interface, the fluid and solid meshes may match node for node, or they may be independent.
A matching mesh makes coupling easier. The same interface nodes carry fluid boundary motion and structural displacement. Kinematic continuity can be imposed directly.
A nonmatching mesh gives more freedom. The fluid and solid can use different element sizes and different refinement patterns. That helps when a thin structure touches a large fluid region. The price is that the interface needs a transfer method, such as:
Lagrange multipliers
Mortar projection
Nitsche’s method
Conservative interpolation
Penalty coupling
The transfer must respect both displacement and force balance. If interpolation loses too much force or maps it inconsistently, the simulation can show artificial energy gain or damping.
Strong and weak interface enforcement
The interface conditions can be imposed strongly or weakly.
Strong enforcement directly ties degrees of freedom together. This works well with matching meshes and compatible unknowns.
Weak enforcement adds the interface conditions to the variational problem. Nitsche’s method is a common example. It can impose constraints without adding multiplier unknowns, and it works well for nonmatching discretizations when the terms are chosen carefully.
Lagrange multipliers introduce extra unknowns that represent interface traction. They are attractive because the coupling force becomes explicit, but the resulting system needs compatible multiplier spaces. Poor choices can produce unstable pressure-like oscillations along the interface.
Choosing finite element spaces changes the behavior
Once the weak form is defined, the next step is selecting finite element spaces for the unknown fields. This choice has a strong effect on accuracy and stability.
The fluid typically needs velocity and pressure spaces that satisfy the incompressibility constraint. The classic warning is the inf-sup, or Ladyzhenskaya-Babuška-Brezzi, condition. In plain terms, not every velocity-pressure pair is stable.
Common choices include:
Fluid element pair | Typical use | Main concern |
Taylor-Hood elements, such as \(P_2/P_1\) | Stable incompressible flow | More degrees of freedom |
Mini elements | Lower-order triangular or tetrahedral meshes | Bubble functions add cost |
Equal-order \(P_1/P_1\) or \(Q_1/Q_1\) with stabilization | Large practical models | Stabilization must be tuned with care |
Discontinuous pressure spaces | Local mass behavior and flexibility | More complex assembly |
The solid side may use displacement-based elements, mixed elements, shell elements, beam elements, or membrane elements. Nearly incompressible solids need attention because low-order displacement elements can lock. Thin structures add another layer because bending behavior is sensitive to element formulation.
For moving fluid meshes, a third field may enter the problem: mesh displacement. The mesh motion is artificial, but it matters. Common choices solve a pseudo-elastic or harmonic extension problem in the fluid mesh:
\[
-\nabla\cdot(k \nabla d_m) = 0
\]
where \(d_m\) is mesh displacement and \(k\) can vary by element size to protect small elements from distortion.
What gets discretized
A typical finite element FSI model may contain these unknowns:
Fluid velocity \(u_f\)
Fluid pressure \(p_f\)
Solid displacement \(d_s\)
Solid velocity \(v_s\), depending on time scheme
Mesh displacement \(d_m\)
Interface multiplier \(\lambda\), if using a multiplier method
After spatial discretization, the continuous problem becomes a set of nonlinear algebraic equations. A simplified monolithic residual may look like
\[
R(U) =
\begin{bmatrix}
R_f(u_f,p_f,d_m) \\
R_c(u_f) \\
R_s(d_s) \\
R_m(d_m,d_s) \\
R_{\Gamma}(u_f,d_s)
\end{bmatrix}
=0
\]
where \(R_f\) represents fluid momentum, \(R_c\) incompressibility, \(R_s\) structural balance, \(R_m\) mesh motion, and \(R_{\Gamma}\) interface coupling.
The exact form changes by method, but the pattern stays the same. FSI discretization creates a coupled system with constraints and moving geometry.
Stabilization is often necessary
Many FSI simulations run at Reynolds numbers where convection matters. Fluid stabilization can reduce spurious oscillations, especially with equal-order velocity-pressure pairs. Common families include SUPG, PSPG, and grad-div stabilization.
The goal is not to hide a bad mesh or an unsuitable time step. Stabilization should improve the discrete problem while keeping the physical behavior credible.
A few practical signs that stabilization or element choice needs review include:
Checkerboard pressure patterns
Artificial interface oscillations
Poor local mass conservation
Sensitivity to small time-step changes
Large residuals near the moving boundary
Solid stabilization can also matter. Thin solids, contact, near-incompressibility, and large rotations can all require special element technology or mixed formulations.
Time discretization and coupling drive stability
Spatial finite elements are only half the story. FSI is time dependent, and the time scheme has a major effect on stability.
For the fluid, common choices include backward Euler, BDF2, and generalized-\(\alpha\) methods. For the structure, Newmark and generalized-\(\alpha\) schemes are common. A good scheme should give enough damping for high-frequency numerical noise without erasing meaningful dynamics.
A simple backward Euler approximation for velocity is
\[
\frac{u^{n+1}-u^n}{\Delta t}
\]
For the structure, displacement and velocity may be updated together, with acceleration included if inertia matters. In strongly coupled FSI, all interface conditions should be evaluated consistently at the same time level.
Monolithic coupling
A monolithic method solves the fluid, solid, mesh motion, and interface equations in one coupled system. Conceptually, this is the cleanest approach.
A block system might have this structure:
\[
\begin{bmatrix}
A_f & B_f^T & C_{fm} & C_{fs} \\
B_f & 0 & 0 & 0 \\
C_{mf} & 0 & A_m & C_{ms} \\
C_{sf} & 0 & 0 & A_s
\end{bmatrix}
\begin{bmatrix}
\Delta u_f \\
\Delta p_f \\
\Delta d_m \\
\Delta d_s
\end{bmatrix}
=
-
\begin{bmatrix}
R_f \\
R_c \\
R_m \\
R_s
\end{bmatrix}
\]
The block names are only schematic. Real systems include nonlinear terms, interface constraints, time terms, and material tangents.
The main strength of monolithic coupling is stability, especially when the fluid and structure densities are close. That case appears in many biomedical and underwater applications. The main cost is implementation complexity. The solver must handle a large, mixed, nonlinear system with very different block behavior.
Newton’s method is often used for the nonlinear solve. Each Newton step needs a linearized system. Accurate tangents improve convergence, but deriving and coding them takes effort.
Partitioned coupling
A partitioned method solves the fluid and solid with separate solvers and exchanges data at the interface. For example:
Use the current structural position to update the fluid mesh.
Solve the fluid problem.
Transfer fluid traction to the structure.
Solve the structural problem.
Transfer displacement or velocity back to the fluid.
Repeat within the time step until the interface is consistent.
This approach is attractive because existing fluid and solid solvers can be reused. It also keeps each physics module easier to manage.
The risk is instability from the added-mass effect. When a light or flexible structure interacts with a dense fluid, simple staggered coupling can fail, even with small time steps. The fluid behaves as if it adds inertia to the structure, and loose coupling may lag that inertia.
Common remedies include:
Subiterations at each time step
Aitken relaxation
Quasi-Newton interface acceleration
Robin interface conditions
Implicit treatment of interface motion
Partitioned methods can perform very well, but the coupling algorithm is part of the numerical method. It cannot be treated as a mere software connection.
Mesh motion is a numerical model of its own
The fluid mesh in ALE FSI must move with the structure on the interface. Inside the fluid domain, its motion is not physical. It is a computational device. Bad mesh motion can ruin an otherwise sound discretization.
A basic harmonic extension spreads boundary displacement into the fluid mesh. It is simple and often works for small motions. Pseudo-elastic mesh motion treats the mesh like a fictitious solid. It can handle larger deformation and allows stiffness to vary across the domain.
For large motion, remeshing may be needed. Remeshing solves element distortion but introduces projection errors because fields must be transferred from the old mesh to the new one. In transient FSI, those transfer errors can affect conservation and interface balance.
Some problems avoid moving meshes by using immersed or fictitious-domain methods. The structure moves through a background fluid mesh, and the interface is captured by constraints or forcing terms. These methods help when motion is very large, but interface accuracy and boundary layer resolution need careful treatment.
The geometric conservation law matters
ALE formulations should respect the geometric conservation law. In simple terms, a moving mesh should not create artificial flow by its motion alone. If the domain moves but the physical solution is uniform, the numerical scheme should preserve that uniformity.
This requirement affects how mesh velocity, time integration, and changing element volumes are computed. Inconsistent geometry updates can introduce small errors that grow over many time steps.
Integration on moving domains
Element integrals are evaluated on the current fluid domain. As the mesh moves, shape-function gradients, Jacobians, and quadrature-point positions change. A code must update these geometric quantities at the right stage of the nonlinear or time iteration.
For nonlinear solids, integration also tracks deformation measures at quadrature points. Hyperelastic materials need stress and tangent evaluations based on the current deformation gradient. Viscoelasticity, plasticity, or damage models add internal variables that must follow the time integration scheme.
A practical discretization workflow
A reliable FSI finite element model usually comes from a sequence of controlled choices. The details vary by application, but the workflow below is a useful baseline.
Define the physics before the mesh
Start with the expected deformation, flow regime, and material behavior.
For the fluid, decide whether incompressible Navier-Stokes is enough or whether compressibility, turbulence modeling, heat transfer, or multiphase effects matter.
For the structure, decide whether the model needs small-strain elasticity, large-deformation hyperelasticity, shells, beams, membranes, or contact.
These decisions determine the weak form and the element types. Mesh density comes later.
Choose compatible spaces
Pick fluid velocity-pressure spaces with stable incompressibility behavior. If using equal-order interpolation, include a tested stabilization method.
Pick structural spaces that match the deformation mode. A thick rubber block, a thin steel panel, and a fabric membrane should not use the same default formulation without review.
For the interface, decide whether the mesh will match. If it will not, choose a conservative and stable transfer method.
Set the time scheme and coupling level
If the structure is heavy and motions are modest, a partitioned staggered method may work well. If the fluid and solid densities are similar or the structure is very flexible, stronger coupling is often needed.
Use time steps based on both fluid and structural scales. A time step that resolves the flow may still miss structural vibration, and the reverse can also happen.
Check conservation and interface balance
Do not judge the model by smooth animations alone. Track numerical measures such as:
Interface displacement mismatch
Interface force balance
Fluid mass conservation
Mechanical energy behavior
Nonlinear residual reduction
Mesh quality over time
These quantities reveal problems early. A visually plausible FSI result can still violate key balances.
Refine with a purpose
Mesh refinement should answer a question. Refine the interface if traction transfer is poor. Refine boundary layers if wall shear matters. Refine the structure where stress gradients are high. Refine in time if phase lag or oscillation frequency changes.
A good convergence study changes one major factor at a time. If space, time, stabilization, and coupling tolerance all change together, the source of improvement stays unclear.
The main takeaway
Finite element FSI succeeds when the discretization treats the interface as the center of the problem. The fluid, structure, mesh motion, and time scheme must agree on how motion and traction pass across that boundary.
A strong model does four things well:
It uses stable finite element spaces for incompressible flow and structural response.
It enforces interface kinematics and dynamics consistently.
It handles mesh motion without corrupting the fluid solve.
It chooses monolithic or partitioned coupling based on the physics, not convenience alone.
The finite element method gives FSI a flexible foundation, but the quality of the result comes from the details: the weak form, the spaces, the coupling, and the checks that prove the numerical exchange is physically credible.
Comments