AFSI: Automated Fluid-Structure Interaction Solver Development for Nonlinear Solid Mechanics

TL;DR

AFSI is an open-source, flexible fluid-structure interaction solver built on FEniCS, capable of simulating large deformations in hyperelastic materials without remeshing, suitable for biomechanical and complex nonlinear applications.

Contribution

It introduces a novel immersed boundary framework within FEniCS for nonlinear FSI, enabling easy customization and efficient simulation of highly deformable structures.

Findings

Successfully simulates large deformations in hyperelastic materials.

Provides a modular Python API for rapid setup and modification.

Achieves efficient performance with hybrid Python/C++ architecture.

Abstract

AFSI is a novel, open-source fluid-structure interaction (FSI) solver that extends the capabilities of the FEniCS finite element library through an immersed boundary (IB) framework. Designed to simulate large deformations in hyperelastic materials (such as cardiac tissue), AFSI avoids the need for expensive remeshing by coupling a Lagrangian representation of the solid with an Eulerian description of the surrounding fluid. This approach retains the full expressiveness of FEniCS's variational formulations, function spaces, and time integration schemes. Implemented in a hybrid Python/C++ architecture, AFSI allows users to define geometries, constitutive models (e.g., the Holzapfel-Ogden law for myocardium), and strain energy functions directly in Python, while delegating performance-critical tasks such as assembly and linear solvers to optimized C++ backends. Its concise and modular Python API facilitates the setup of FSI simulations, enabling users to easily modify discretization strategies or analyze results using standard FEniCS post-processing tools. By combining the flexibility of FEniCS with a robust immersed boundary formulation, AFSI empowers rapid prototyping of complex nonlinear solid-fluid interaction problems, making it a powerful tool for simulating biomechanical systems and other applications involving highly deformable structures in flow.

Tables

Table 1. Table 1 : Parameters for the idealized heart valve model with isotropic properties.

Symbol	Value	Unit
$T$	3.0	s
$\rho$	100	$\frac{\text{g}}{{\text{cm}}^3}$
$\mu$	10.0	$\text{g}/\left(\text{cm}\cdot \text{s}\right)$
$E$	$5.6\times {10}^{-7}$	$\frac{\text{dyn}}{{\text{cm}}^2}$
$\nu$	0.4
$\Delta t$	$1\times {10}^{-4}$	s

Table 2. Table 2 : Parameters for the idealized heart valve model with anisotropic properties.

Symbol	Value	Unit
$T$	3.0	s
$\rho$	100	$\frac{\text{g}}{{\text{cm}}^3}$
$\mu$	10.0	$\text{g}/\left(\text{cm}\cdot \text{s}\right)$
$C_0$	$2.0\times {10}^7$	Pa
$C_1$	$5\times C_0$	Pa
${\kappa }_s$	$1.0\times {10}^7$	Pa
$\Delta t$	$1\times {10}^{-4}$	s

Table 3. Table 3 : Parameters for the idealized left ventricle model.

Symbol	Value	Unit
$T$	3.0	s
$\rho$	1	$\frac{\text{g}}{{\text{cm}}^3}$
$\mu$	0.04	$\text{g}/\left(\text{cm}\cdot \text{s}\right)$
$C$	$10$	kPa
$b_f$	$1$	-
$b_t$	$1$	-
$b_{fs}$	$1$	-
$\Delta t$	$1\times {10}^{-4}$	s

Equations

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}(\mathbf{x},t)+\mathbf{u}(\mathbf{x},t)\cdot \nabla \mathbf{u}(\mathbf{x},t)\right)$$

$$\nabla \cdot \mathbf{u}(\mathbf{x},t)$$

$${\int }_{B_r}\mathbf{F}(\mathbf{X},t)\cdot \mathbf{V}(\mathbf{X})d\mathbf{X}$$

$$\mathbf{f}(\mathbf{x},t)$$

$$\frac{\partial \mathcal{X}(\mathbf{X},t)}{\partial t}$$

$$\mathbf{u}(\mathbf{x},t)$$

$$\mu \frac{\partial \mathbf{u}(\mathbf{x},t)}{\partial \mathbf{n}}-p(\mathbf{x},t)\mathbf{n}$$

$$\mathbf{u}(\mathbf{x},0)$$

$$\mathbb{P}(\mathbf{X},t)=\frac{\partial \Psi (\mathbb{F}(\mathbf{X},t))}{\partial \mathbb{F}(\mathbf{X},t)},$$

$$\mathbb{F}(\mathbf{X},t)=\frac{\partial \mathcal{X}(\mathbf{X},t)}{\partial \mathbf{X}},$$

$$\Psi (\mathbb{F})=W_{\text{iso}}(\mathbb{F})+U_{\text{vol}}(J),$$

$$W=\frac{1}{2}{\mu }_s(I_1-d)-{\mu }_s\ln (J)+\frac{1}{2}\lambda \ln (J{)}^2,$$

$$\mathbf{v}(\mathbf{x},t)=\{\begin{array}{ll}5(\sin (2\pi t)+1.1)y(1.61-y){\mathbf{e}}_x,& t>0,\\ 0,& \text{ otherwise },\end{array}$$

$$\Psi =\frac{{\mu }^s}{2}\left(\text{tr}(\mathbb{C}-3)\right)-{\mu }^s\text{ln}(J)+\frac{{\lambda }^s}{2}{\left(\text{ln}(J)\right)}^2,$$

$$\Psi =\frac{C_0}{2}\left({\bar{I}}_1-3\right)+C_1\left(\text{exp}({\bar{I}}_{4f}-1)-{\bar{I}}_{4f}\right)+\frac{{\kappa }_s}{2}\left(\frac{1}{2}\left(J^2-1\right)-\text{ln}\left(J\right)\right)$$

$$\mathbf{X}(u,v)=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{c}{r}_s\sin u\cos v\\ r_s\sin u\sin v\\ r_l\cos u\end{array}\right).$$

$$\Psi =\frac{C}{2}\left(\text{exp}(Q)-1\right)$$

$$Q=b_f{E}_{11}^2+b_t(E_{22}^2+E_{33}^2+E_{23}^2+E_{32}^2)+b_{fs}(E_{12}^2+E_{21}^2+E_{13}^2+E_{31}^2),$$

Full text

Keywords: Fluid-structure interaction, immersed boundary finite element method, automated finite element modeling, nonlinear solid mechanics

1 Introduction

Simulating the interactions between flexible solid materials and surrounding fluids is of fundamental importance across a wide range of engineering and biomedical applications. For instance, evaluating the blood bumping function of ventricle, the open and close of heart valves during cardiac cycles [1, 2, 3, 4], and the motion of soft robotic swimmers [5] or biomimetic underwater robots [6, 7, 8] all require a detailed understanding of fluid-structure interaction (FSI). These systems typically involve large deformations, nonlinear material responses, and intricate fluid dynamics, posing significant challenges for both experimental measurement and numerical modeling [9, 10].

Physical experiments in FSI domains are often costly and constrained by limited spatial and temporal resolution. As a result, accurate and efficient numerical FSI solvers have become indispensable tools for understanding and predicting the behavior of these complex coupled systems. However, despite the availability of numerous scientific computing frameworks, practical FSI simulation still faces significant challenges. The classical Arbitrary Lagrangian–Eulerian (ALE) method is known for its numerical stability and accuracy [11], but it often requires generating conforming fluid meshes that move and deform with the solid geometry, which is an operation that becomes prohibitively expensive in three-dimensional simulations. As observed in the simulation of prosthetic aortic valves, the remeshing overhead in ALE-based FSI can far exceed the cost of solving the coupled system itself, severely limiting usability in dynamic scenarios [12]. In contrast, the immersed boundary (IB) method employs nonconforming meshes between fluid and structure. Although its application to materials with finite volume and nonlinear material responses developed later than ALE methods [13], the IB method has rapidly advanced, particularly in computational biology, and is now widely used to simulate biological fluid dynamics, especially in cardiac mechanics [14]. It supports nonconforming discretizations of fluid and structure, eliminating the need for dynamically generated body-fitted meshes. This is especially advantageous in problems involving large deformations, structural displacements, or contact interactions [15]. Although many IB-based FSI solvers exist, only a few support a wide range of hyperelastic material models. IBAMR is one of the most sophisticated IB-based solvers. However, its complexity and heavy reliance on external libraries such as SAMRAI, PETSc, and libMesh—combined with its extensive low-level C++ codebase—pose a steep learning curve for new users.

In our work, we address these challenges through a combination of tailored numerical methods and efficient software design. We adopt a nodal immersed boundary finite element (IBFE) method [16], which decomposes the FSI problem into an Eulerian Navier-Stokes solver for the fluid and a Lagrangian solid solver. The coupling between these two domains is achieved via integral transforms that transfer information between Eulerian and Lagrangian variables. A key contribution of our framework is the C++ implementation of these coupling operators, which uses shared-memory parallelism via Intel TBB, and is exposed to Python through nanobind. This design allows users to write simulations in compact, human-readable Python code, typically within a few hundred lines per example, where each line closely mirrors the corresponding mathematical expression. Users do not need to manage or even be aware of the underlying parallelism, which is fully abstracted away. Besides, it is integrated into FEniCS, which translates high-level variational forms into optimized C++ code at runtime. This eliminates the need for manual low-level finite element coding and significantly streamlines the development process. FEniCS also provides built-in support for parallel computing, enabling efficient large-scale simulations. The framework is designed to be modular and extensible, allowing flexible experimentation with different finite element discretizations, time-stepping schemes, and nonlinear material models. Its effectiveness is demonstrated through several numerical experiments, including blood flow-driven motion of a two-dimensional valve and torsion of an elastic beam immersed in fluid. These examples showcase the solver’s ability to handle complex geometries, large deformations, and nonlinear constitutive behaviors.

The remainder of this paper is structured as follows. Section 2 presents the governing equations of the coupled fluid-structure system, hybrided with Eulerian and Lagrangian frameworks, constitutive models for hyperelastic materials, and the nodal IBFE method formulation. Section 3 describes the implementation of the solver, with a focus on temporal discretization, the automated finite element framework, and its hybrid Python/C++ architecture. Section 4 demonstrates the solver’s accuracy and robustness through benchmark cases. Finally, Section 5 summarizes the main contribution of this study and discusses future research directions.

2 Methodology

2.1 Fluid-structure interaction problem

This study is based on the IB method derived by Boffi et al.[17], which is used to solve FSI systems. The method is mathematically described by the following system:

[TABLE]

The variables and components in the above system are defined as follows:

• •

$\rho$ is the constant density shared by the fluid and structure.

• •

The velocity $\mathbf{u}(\mathbf{x},t)$ and pressure $p(\mathbf{x},t)$ are defined on the fixed Eulerian computational domain $\Omega$ and are governed by Eq. (1a)–Eq. (1b).

• •

$\mu$ is the dynamic viscosity of the incompressible Newtonian fluid.

• •

The external force density $\mathbf{f}(\mathbf{x},t)$ acts as a coupling term between the Lagrangian and Eulerian frameworks (Eq. (1d)).

• •

The Lagrangian elastic force density $\mathbf{F}(\mathbf{X},t)$, defined on the reference configuration $B_{r}$ of the immersed structure, satisfies the weak formulation in Eq. (1c) for all smooth test functions $\mathbf{V}(\mathbf{X})$. Here, $\mathbb{P}(\mathbf{X},t)$ is the first Piola–Kirchhoff stress tensor.

• •

The structure’s motion is described by its current position $\mathcal{X}(\mathbf{X},t)$ and velocity $\mathbf{U}(\mathbf{X},t)$, which are related to the surrounding fluid velocity via interpolation using the Dirac delta function (Eq. (1e)).

• •

The Eulerian spatial coordinates are denoted by $\mathbf{x}=(x_{1},\dots,x_{d})\in\Omega$, and the Lagrangian coordinates by $\mathbf{X}=(X_{1},\dots,X_{d})\in B_{r}$.

• •

The Dirac delta function $\delta(\mathbf{x}-\mathcal{X}(\mathbf{X},t))$ serves to mediate the two-way coupling between fluid and structure by enabling force spreading and velocity interpolation between the Eulerian and Lagrangian frames.

The computational model is subject to the following initial and boundary conditions:

[TABLE]

where $\partial\Omega=\partial\Omega_{D}\cup\partial\Omega_{N}$ is the boundary of the fluid domain, $\mathbf{w}_{D}(\mathbf{x},t)$ is the prescribed velocity on the boundary, $\mathbf{w}_{N}(\mathbf{x},t)$ is the prescribed traction on the boundary, and $\mathbf{u}_{0}(\mathbf{x})$ is the initial velocity field.

2.2 Mechanical properties of the immersed structure

This paper focuses on hyperelastic structures, whose material response is described by a strain energy density function $\Psi(\mathbb{F}(\mathbf{X},t))$. Accordingly, the first Piola–Kirchhoff stress tensor $\mathbb{P}(\mathbf{X},t)$ is computed as

[TABLE]

where $\mathbb{F}(\mathbf{X},t)$ is the deformation gradient, defined by

[TABLE]

and its determinant is denoted by $J(\mathbf{X},t)=\det\left(\mathbb{F}(\mathbf{X},t)\right)$.

Since the solid materials considered in this study are incompressible or nearly incompressible, we employ a dilatational penalty method to model their mechanical behavior (see [18, 19]). In this approach, the strain energy density function $\Psi(\mathbb{F})$ is decomposed into two parts: an isochoric component and a volumetric component:

[TABLE]

representing the volume-preserving and volume-changing elastic response of the material, respectively.

3 Implementation

3.1 Discretization

For the system Eq. (1), the time interval $[0,T]$ is uniformly divided into $N$ non-overlapping subintervals $(t_{n},t_{n+1}]$, where $n=1,2,\dots,N$, and the time step size is defined as $\Delta t=t_{n}-t_{n-1}$. The notation $(\cdot)^{n}$ denotes the value of a variable at time step $t_{n}$. some varibales use the n+1, some variables use n.

The background domain is triangulated using quadrilateral elements for two-dimensional problems and hexahedral elements for three-dimensional problems. The velocity field is approximated using $Q2$ finite elements, while the pressure field is approximated using $Q1$ elements. The interpolation nodes for the finite element basis functions are illustrated in Figures 1 and 2.

The structural domain is triangulated using triangular elements for two-dimensional problems and tetrahedral elements for three-dimensional problems. The structural force field $\mathbf{F}(\mathbf{X},t)$ and current position field $\mathcal{X}(\mathbf{X},t)$ are approximated using $P2$ finite elements. The interpolation nodes for the finite element basis functions are illustrated in Figures 3 and 4.

3.2 Automated finite element solution

The FEniCS framework provides a convenient way to express strain energy density functions in symbolic form, maintaining a close correspondence with their mathematical notation. For example, the strain energy density function of a neo-Hookean material is given by

[TABLE]

where the notations have been defined above. This expression can be implemented in Python within FEniCS as follows.

⬇

1def NeoHookeanModel(F, mu_s, lmbda, mesh):

2 d = mesh.geometry().dim()

3 C = F.T*F

4 J = det(F)

5 I1 = tr(C)

6 W = (mu_s/2)(I1 - d) - mu_sln(J) + (lmbda/2)*(ln(J))**2

7 return W

From this symbolic form, the first Piola-Kirchhoff stress tensor in the weak formulation Eq. (1c) can be derived automatically in FEniCS simply by calling "diff(W, F)". As the discretization method defined in subsection 3.1, FEniCS generates and compiles highly optimized C++ code at runtime, ensuring both clarity in implementation and efficiency in computation.

As for the solver for the incompressible Navier-Stokes equations, Eq. (1a) and Eq. (1b), the discretization approach has already been described in Subsection 3.1. As discussed in [20], various implementations are possible; in our work, unless stated otherwise, we adopt Chorin’s projection method.

3.3 Eulerian-Lagrangian coupling

The Eulerian-Lagrangian coupling is implemented through the force spreading (Eq. (1d)) and velocity interpolation (Eq. (1e)) operators. While these formulations are not natively supported in FEniCS due to their requirement for interaction between distinct finite element spaces, we have developed a novel framework that integrates this capability directly within FEniCS, enabling complete solution of the immersed boundary system Eqs. (1). This integration represents a key contribution of our work.

Our implementation builds upon the nodal coupling scheme introduced by [16], with adaptations for finite element discretization. In [16], the scheme is implemented using a finite difference method with Eulerian interaction nodes defined on a staggered grid. In the present work, the coordinates of the degrees of freedom illustrated in Figures 1 and 2 chosen as the Eulerian interaction nodes. The Lagrangian interaction nodes follow the arrangement in Figures 3 and 4, as in [16].

The computational core is implemented in C++ for optimal performance, with Python interfaces exposed via nanobind [21]. This hybrid approach preserves computational efficiency while maintaining compatibility with FEniCS’s Python ecosystem. It enables seamless integration with FEniCS’s existing infrastructure and user-friendly Python interface. The implementation provides two key classes:

• •

IBMesh: Creates structured background meshes and manages DoF mappings for efficient coupling operations;

• •

IBInterpolation: Handles the core coupling operations between Eulerian and Lagrangian domains.

They can be initialized using the following code.

⬇

1eulerian_mesh = IBMesh([Point(0, 0), Point(1, 1)], [10, 10], order_velocity)

2lagrangian_mesh = RectangleMesh(Point(0.2, 0.2), Point(0.8, 0.8), 10, 10)

3ib = IBInterpolation(eulerian_mesh, lagrangian_mesh)

The key coupling operations of Eq. (1d) and Eq. (1e) are implemented through the following interface.

⬇

1ib.fluid_to_solid(flow_velocity._cpp_object, solid_velocity._cpp_object)

2ib.solid_to_fluid(flow_force._cpp_object, solid_force._cpp_object)

3.4 Parallelism

[22] Gathered on the same rank

currently the bottleneck of the entire problem

4 Examples

To demonstrate the applicability and reliability of the developed software AFSI in addressing FSI problems in biomechanics, two representative numerical benchmark cases are examined in this section. The first is a two-dimensional heart valve model, designed to evaluate the software’s capability in handling strongly coupled interactions between thin structures and surrounding fluid. The second is a three-dimensional idealized left ventricle model, aimed at assessing the stability and accuracy of AFSI under complex geometries and nonlinear coupling conditions. It is important to note that these benchmark cases are intended solely for research validation and will not be included in the public release of the software. Instead, a set of simplified and reproducible example cases will be provided at release, serving as references for users to understand the core functionalities and basic usage of AFSI.

4.1 Idealized heart valve

Here, a classical benchmark problem of an idealized two-dimensional heart valve is adopted, which has been extensively investigated and discussed in the literature [23, 24, 25, 26]. Based on the setup proposed in [23, 24], a thin cantilever beam—representing a valve leaflet—was fixed to each of the upper and lower walls of a channel filled with an incompressible Newtonian fluid, as shown in Figure 1. This configuration was used to assess the accuracy of AFSI in capturing heart valve dynamics under pulsatile flow conditions. To avoid structural contact, the beams were intentionally designed to be slightly shorter than the distance between the channel centerlines. This setup can also be interpreted as a simplified model of a check valve.

For this problem, the boundary conditions are set as follows: The inlet is located on the left boundary of the domain $\Omega$, where the fluid is driven by a time-dependent velocity profile

[TABLE]

in order to replicate the characteristics of pulsatile flow. No-slip boundary conditions are imposed on the upper and lower walls of the domain, corresponding to zero velocity ($\mathbf{u}=\mathbf{0}$). The right boundary serves as the outlet and is assigned a free traction condition to allow the fluid to exit the domain freely.

Under the same geometry and boundary conditions, we next analyze two valve materials with different constitutive models: isotropic and anisotropic.

4.1.1 FSI validation: leaflets with isotropic constitutive Materials

The first case considers the leaflets as isotropic hyperelastic materials, modeled using the Neo-Hookean model with a dilatational penalty term, and the model is defined as follows:

[TABLE]

where $\mu^{s}=\frac{E}{2(1+\nu)}$ and $\lambda^{s}=\frac{E\nu}{(1+\nu)(1-2\nu)}$ are the Lamé parameters. The parameter settings used for the simulation of this model are listed in Table 1. It should be noted that the constitutive model originally represents a compressible material; however, due to the chosen Poisson’s ratio of $\nu=0.4$, the material is considered nearly incompressible in practice.

In this case, the structural domain is discretized into $4\times 40$ quadrilateral elements, while the fluid domain is discretized using a structured quadrilateral mesh of size $256\times 64$. Figure 2 shows the velocity component in the x-direction and the displacement magnitude of the valve leaflets at different times during the cardiac cycle from 2 s to 3 s. The results exhibit excellent agreement with the simulations reported in [23]. As illustrated in the figure, when the flow rate approaches its peak, a high-velocity jet forms on the upstream surface of the valve, while a flow stagnation region develops immediately downstream of the leaflets. Figure 3 presents a comparison of the time-dependent displacements at the tip of the upper valve leaflet in the $x$ and $y$ directions. The results from this study are compared with those obtained using the ALE-FSI method [24] and the distributed Lagrange multiplier (DLM) type FSI method reported in [23]. It is worth noting that the ALE-FSI method employs the Saint-Venant–Kirchhoff constitutive model, whereas the DLM method adopts the same material model as used in the present study. As shown in the figure, the displacement curves produced by the proposed method exhibit excellent agreement with both reference solutions, thereby validating the accuracy of the current approach.

4.1.2 FSI validation: leaflets with anisotropic constitutive Materials

In reality, heart valves are biologically complex and compositionally heterogeneous tissues that exhibit pronounced anisotropic mechanical behavior. To more accurately capture their biomechanical characteristics, the second case models the valve leaflets as anisotropic hyperelastic materials, employing the anisotropic constitutive model proposed in [23], whose specific formulation is given below,

[TABLE]

where $C_{0}$ and $C_{1}$ are the material parameters. The first term in Eq. (3) captures the isotropic response of the matrix, whereas the second term accounts for a nonlinear stiffening effect along a preferred direction, resulting in an overall transversely isotropic mechanical behavior. The parameter settings used for The anisotropic material model are listed in Table 2.

In this case, the structural and fluid domains adopt the same mesh configuration as in the first case. Three different fiber orientation angles ($45^{\circ}$, $60^{\circ}$, and $75^{\circ}$), as illustrated in Figure 4, are considered to assess the capability of the AFSI software in modeling anisotropic constitutive materials. The simulation results indicate that the chosen material properties and constitutive model produce a more compliant isotropic background matrix in the leaflets, leading to more complex deformation at certain fiber orientations. Figure 5 shows the displacements of the upper leaflet tip in the $x$ and $y$ directions during the final cardiac cycle. Figure 6 illustrates the motion of the valve leaflets with different fiber orientations over the same period. As shown, the displacement and overall dynamic behavior of the leaflet with a $60^{\circ}$ fiber orientation lie between those observed at $45^{\circ}$ and $75^{\circ}$. Among the three cases, the $45^{\circ}$ orientation yields the largest leaflet tip displacement, while the $75^{\circ}$ orientation results in the smallest. These results underscore the significant role of fiber orientation in influencing the mechanical response of the valve leaflets, which is crucial for evaluating valve performance under different physiological conditions.

4.2 Idealized left ventricle

To demonstrate the numerical simulation capabilities of the AFSI software in handling complex geometries and structural models, a three-dimensional idealized left ventricle is employed as a test case. This model serves as a simplified representation of the human left ventricle and is widely used as a benchmark for verifying the accuracy of cardiac mechanics solvers [27, 28]. Figure 7 depicts the initial geometry of the idealized left ventricle, along with the associated fixed computational domain. The undeformed geometry is defined through the following parameterization of a truncated ellipsoid:

[TABLE]

The boundary parameters are defined as follows:

• •

Endocardial surface: $r_{s}=0.7\leavevmode\nobreak\ \text{cm}$, $r_{l}=1.7\leavevmode\nobreak\ \text{cm}$, $u\in\left[-\pi,-\arccos\frac{5}{17}\right]$, and $v\in[-\pi,\pi]$;

• •

Epicardial surface: $r_{s}=1\leavevmode\nobreak\ \text{cm}$, $r_{l}=2.0\leavevmode\nobreak\ \text{cm}$, $u\in[-\pi,-arccos\frac{5}{20}]$, and $v\in[-\pi,\pi]$;

• •

Base plane: $z=0.5\leavevmode\nobreak\ \text{cm}$, which is implicitly defined by the given range of $u$.

The passive mechanical behavior of the left ventricle is modeled using a transversely isotropic hyperelastic constitutive law proposed by Guccione et al. [29]:

[TABLE]

where

[TABLE]

and $\mathbb{E}=\frac{1}{2}(\mathbb{F}^{T}\mathbb{F}-\mathbb{I})=(E_{ij})$ is the Green-Lagrange strain tensor. Here, $C$, $b_{f}$, $b_{t}$, and $b_{fs}$ are material parameters. The geometry of the idealized left ventricle is discretized using 17,879 tetrahedral elements, and the fluid computational domain $\Omega$ is resolved on a $64\times 64\times 64$ Cartesian mesh. The parameters employed in the simulation are summarized in Table 3. A uniform pressure load of $10\leavevmode\nobreak\ \text{kPa}$ is applied to the endocardial surface to simulate physiological loading conditions.

5 Conclusions

This implementation enables seamless integration of the IB method within the FEniCS framework for FSI simulations, offering several notable advantages. First, the nodal coupling scheme naturally accommodates arbitrary high-order finite element spaces, ensuring both flexibility in spatial discretization and high-order accuracy. Second, the methodology completely avoids the need for remeshing when handling large deformations and displacements, thereby substantially reducing preprocessing complexity and computational overhead. Third, the framework permits straightforward selection among various time discretization schemes and solution algorithms through FEniCS’s native time-stepping capabilities. Finally, the implementation supports automatic derivation and implementation of finite element formulations directly from strain energy functions, which particularly benefits modeling of solids with complex constitutive behavior and significantly streamlines the development of new material models.

Author contributions

Pengfei Ma: Methodology (lead); Software(lead); Writing - original draft (lead). Li Cai: Supervision (lead); Writing - review & editing (equal). Xuan Wang: Visualization (lead); Writing - original draft (supporting); Writing - review & editing (equal). Hao Gao: Writing - review & editing (equal).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

Pengfei Ma, Li Cai and Xuan Wang were supported by the National Natural Science Foundation of China (Grant No. 12271440). Hao Gao was supported by the British Heart Foundation (PG/22/10930) and the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom (EP/S030875/1).