Finite Element Method for Strongly Coupled Fluid-Structure Interaction of a Vertical Flap in a Channel and Aneurysm Hemodynamics

1. Introduction

Fluid-structure interaction (FSI) is a phenomenon that occurs when a fluid interacts with a solid object. This interaction can be studied by examining changes in the structural properties of the solid and the complexity of the fluids. The key parameters that play a significant role in FSI include fluid velocity, fluid viscosity, Young’s modulus of the solid, density of the solid, and the Poisson ratio of the structure [1,2]. The study of FSI and different multi-physics problems has received increasing attention in recent years due to their importance in various fields of applied sciences and engineering, including turbo machinery [3], bridge decks [4], aerospace [5,6,7], parachute dynamics [8], flow around elastic structures, blood flow in the cardiovascular system [9], the study of hemodynamics [10,11,12], and the flutter analysis of airplanes [13].

Recent developments in this field include the shift from defined purpose to more general approaches, the need to capture complex and general systems, and the demand for robust high-quality approaches that can accurately predict outcomes in complex cases. As a result, there has been a growing interest in the development and application of simulation techniques for FSI in recent decades [14,15,16].

The FSI interface utilizes an Arbitrary Lagrange Euler method, where the fluid flow is formulated using an Eulerian description in a spatial frame, while solid mechanics is formulated using a Lagrangian description, i.e., a material frame. Several researchers have proposed different techniques to address the FSI problem, including an extended ALE method proposed by S. Basting et al. [17] and various mesh moving techniques for monolithically-coupled fluid-structure interactions in an Arbitrary Lagrangian-Eulerian coordinate analyzed by T. Wick [18]. J. Hron, M. Razzaq, and S. Turek analyzed a monolithic algorithm to solve the problem of time-dependent interaction between an incompressible viscous fluid interacting with an elastic solid, which is the method of choice [1,19].

This study focuses on the deformations of an elastic vertical flap in a fluid channel for various inlet mean velocities. Analysis of the displacements of the elastic wall geometries and pressure distribution is conducted, and changes in the flow patterns are observed with the deformations. Additionally, the proposed model is applied to hemodynamical geometries and wall shear stress analysis to demonstrate the significant benefits of this method. The results of this study have broad implications, including great health importance and extensive mathematical model incorporation. Furthermore, the FSI case presented in this study has applications in MEMS devices, where the deformation of the flap due to fluid interaction can provide useful insights about structural properties and aid in designing transduction mechanisms for these devices [20].

The remainder of this article is organized as follows. Section 2 describes the mathematical formulation of the FSI problem. Section 3 explains the FEM discretization of the problem, and Section 4 outlines the solution algorithm procedures used to numerically solve the problem. The simulated results for the flap geometry and aneurysm configuration are presented in Section 5, followed by conclusions and future outlooks and suggestions in Section 6.

2. FSI Problem Formulation

A universal FSI problem comprises of fluid description, solid description, interface conditions and constraints for remaining boundaries. In this article, the incompressible Newtonian fluid flow interacting with a vertical elastic flap and elastic boundary are taken.

2.1. Fluid Model

The fluid properties are represented by the velocity ${\mathbf{v}}^f$ and pressure$p^f$, where the fluid is assumed to be Newtonian with constant density, ${\rho }^f$, and kinematic viscosity, ${\nu }^f$. The balance equation is given by:

$$\begin{array}{cc}\hfill {\varrho }^f\frac{\partial v^f}{\partial t}+{\varrho }^f(\nabla v^f)(v^f-\frac{\partial u^f}{\partial t})& =\mathrm{div}{\mathbf{o}\mathbf{e}}^f\mathrm{in}{\Omega }_t^f,\hfill \\ \hfill \mathrm{div}{v}^f& =0\mathrm{in}{\Omega }_t^f\hfill \end{array}.$$

(1)

In order to solve the balance equation, it is crucial to establish the fundamental relationship for the stress tensor. A constant density Newtonian fluid is used for this purpose, which is expressed as:

$$\begin{array}{c}\hfill {\mathbf{o}\mathbf{e}}^f=-p^f\mathbf{I}+\mu (\nabla {\mathbf{v}}^f+{(\nabla {\mathbf{v}}^f)}^T),\end{array}$$

(2)

here, $\mu$ represents the dynamic viscosity and $\mathbf{I}$ denotes the identity matrix. The Lagrangian multiplier, denoted by $p^f$, corresponds to the incompressibility constraint in (1). The stress tensor is represented by ${\mathbf{o}\mathbf{e}}^f$. The negative term $-p^f\mathbf{I}$ refers to the inviscid reactive component of the Cauchy stress tensor.

2.2. Structure Model

Let us consider a solid object that is elastic and has variable density. Let ${\mathbf{u}}^{\mathbf{s}}$ be the displacement and ${\mathbf{v}}^{\mathbf{s}}$ be its velocity. The balance equation can be expressed as follows:

$${\rho }^s\frac{\partial {\mathbf{v}}^s}{\partial t}+{\rho }^s(\nabla {\mathbf{v}}^s){\mathbf{v}}^s=\mathrm{div}{\sigma }^s,\mathrm{in}{\Omega }_t^s,$$

(3)

where ${\sigma }^s$ denotes the stress tensor. Alternatively, we can write the above equation in a Lagrangian description with respect to a certain initial state ${\Omega }^s$:

$${\rho }^s\frac{{\partial }^2{\mathbf{u}}^{\mathbf{s}}}{\partial t^2}=\mathrm{div}\left(J{\sigma }^s{F}^{-T}\right),\mathrm{in}{\Omega }^s,$$

(4)

here, $F^{-T}={\left(F^{-1}\right)}^T$, $F^{-1}$ is the inverse deformation gradient, and $J=detF$. The equations for the stress tensor for compressible and incompressible structures can be treated in the same way. The undeformed structure density is denoted by ${\rho }^s$. Poisson ratio ${\nu }^s$ and Young’s modulus E determine the elasticity of the structure as follows:

$${\nu }^s=\frac{{\lambda }^s}{2\left({{\mu }^s+\lambda }^s\right)}E=\frac{{\mu }^s(2{\mu }^2+3{\lambda }^s)}{({\lambda }^s+{\mu }^s)},$$

(5)

where ${\mu }^s$ and ${\lambda }^s$ are defined as:

$${\mu }^s=\frac{E}{2({\nu }^s+1)}{\lambda }^s=\frac{E{\nu }^s}{({\nu }^s+1)(1-2{\nu }^s)}.$$

(6)

The Poisson ratio for incompressible structures is $1/2$, while for compressible structures, it is less than $1/2$. The relationship between stress and strain relies on the second Piola-Kirchhoff stress tensor S and the Green Lagrangian strain tensor E, which are functions of the Green Lagrangian strain tensor. The second Piola-Kirchhoff stress tensor, as described in [2], can be obtained from the Cauchy stress tensor ${\sigma }^s$ as follows:

$$S^s=J{F}^{-1}{\sigma }^s{F}^{-T},$$

(7)

where J is the determinant of the deformation gradient tensor F (see [21]), given by

$$F=I+\nabla {\mathbf{u}}^{\mathbf{s}}.$$

(8)

The Green Lagrangian strain tensor E can be expressed as

$$E=\frac{1}{2}(F^{T}F-I),$$

(9)

the St. Venant-Kirchhoff material model defines the Cauchy stress tensor ${\sigma }^s$ in this current study as:

$${\sigma }^s=\frac{1}{J}\mathbf{F}({\lambda }^s\left(trE\right)I+2{\mu }^{s}E){\mathbf{F}}^T,{\mathbf{S}}^s={\lambda }^s\left(trE\right)I+2{\mu }^{s}E,$$

(10)

where J denotes the determinant of the deformation gradient tensor F (see [21]), which is expressed as:

$$F=I+\nabla {\mathbf{u}}^{\mathbf{s}}.$$

(11)

Mooney Rivlin Modeling Using Polynomial-Based Functions

A polynomial-based hyperelastic model has been implemented using the COMSOL multiphysics module [22], with two material model parameters. The model’s equation for strain energy is defined as follows:

$$W=\sum _{i,j=0}^n{C}_{i,j}{(I_1-3)}^i{(I_2-3)}^j+{\displaystyle \frac{1}{2}}K{(J_{el}-1)}^2.$$

(12)

The strain energy density definition for the Mooney Rivlin model with two parameters is given by the following two equations [10]:

$$W_{siso}=C_{1,0}(I_1-3)+C_{0,1}(I_2-3),W_{svol}={\displaystyle \frac{1}{2}}K{(J_{el}-1)}^2.$$

(13)

The material properties used for the simulation are the bulk modulus of $5.6Pa$, and coefficients of the model $C_{1,0}=75Pa$ and $C_{0,1}=100Pa$. The fluid-structure interaction problem is constructed using a pseudo-solid mapping method [23]. The resulting dimensionless system, incorporating the material relation previously mentioned, is presented below:

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{u}}{\partial t}=& \left\{\begin{array}{cc}\Delta \mathbf{u}\hfill & \mathrm{in}{\Omega }^f,\hfill \\ \mathbf{v}\hfill & \mathrm{in}{\Omega }^s,\hfill \end{array}\right.\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{v}}{\partial t}=& \left\{\begin{array}{c}-(\mathbf{v}-\frac{\partial \mathbf{u}}{\partial t}){\mathbf{F}}^{-1}\left(\mathrm{Grad}\mathbf{v}\right)\hfill \\ +\mathrm{Div}\left(-J{p}^f{\mathbf{F}}^{-T}+J\mu \mathrm{Grad}\mathbf{v}{\mathbf{F}}^{-1}{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\Omega }^f,\hfill \\ \frac{1}{\beta }\mathrm{Div}\left(-J{p}^s{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\Omega }^s,\hfill \end{array}\right.\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill 0=& \left\{\begin{array}{cc}\mathrm{Div}\left(J\mathbf{v}{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\Omega }^f,\hfill \\ J-1\hfill & \mathrm{in}{\Omega }^s,\hfill \end{array}\right.\hfill \end{array}$$

(16)

the parameter $\beta$ represents the ratio of solid to fluid density, i.e., $\beta =\frac{{\varrho }^s}{{\varrho }^f}$. In order to maintain the no-slip conditions at the fluid structure interface, the following Neumann and Dirichlet constraints are applied.

$$\begin{array}{cc}\hfill {\mathbf{o}\mathbf{e}}^f\mathbf{n}& ={\mathbf{o}\mathbf{e}}^s\mathbf{n}\hfill \\ \hfill {\mathbf{v}}^f& ={\mathbf{v}}^s.\hfill \end{array}$$

(17)

3. FEM Discretization

The Taylor-Hood pair $P_2{P}_1$ is a popular finite element space selection due to its optimal convergence properties and satisfaction of the LBB stability condition [24]. By transitioning from a discontinuous space to a continuous one, the pressure space’s dimension is reduced, thereby saving degrees of freedom. Figure 1 illustrates the degrees of freedom.

The basis function ${\varphi }_i$ of $P_1$ element is given by [25]

$${\varphi }_i\left(x\right)=a_i+b_{i}x+c_{i}y$$

(18)

and the basis vector set is

$$\langle 1,x,y\rangle .$$

(19)

For $P_2$ element, we have

$${\varphi }_i\left(x\right)=a_i+b_{i}x+c_{i}y+d_i{x}^2+e_i{y}^2+f_{i}xy$$

(20)

and the basis vector set is

$$\langle 1,x,y,x^2,xy,y^2\rangle$$

(21)

Weak Formulation

To discretize in space over the time interval $I=[0,T]$, the standard Galerkin finite element method is applied. The equations (14)–() are multiplied by the test functions $\zeta ,\xi ,\gamma$, and integrated over the space and time interval I. By applying integration by parts on some of the terms and incorporating the boundary conditions, the following outcomes are obtained:

(22)

(23)

$$\begin{array}{cc}\hfill 0=& {\int }_0^T{\int }_{{\Omega }^s}(J-1)\gamma dVdt+{\int }_0^T{\int }_{{\Omega }^f}\mathrm{Div}\left(J\mathbf{v}{\mathbf{F}}^{-T}\right)\gamma dVdt.\hfill \end{array}$$

(24)

Upon applying FEM space discretization, we obtained a nonlinear algebraic system of equations, which can be expressed as follows:

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{\mathbf{S}}_{\mathbf{u}\mathbf{u}}& {\mathbf{S}}_{\mathbf{u}\mathbf{v}}& \mathbf{0}\\ {\mathbf{S}}_{\mathbf{v}\mathbf{u}}& {\mathbf{S}}_{\mathbf{v}\mathbf{v}}& k\mathbf{B}\\ c_{\mathbf{u}}{\mathbf{B}}_{\mathbf{s}}^{\mathbf{T}}& c_{\mathbf{v}}{\mathbf{B}}_f^T& \mathbf{0}\end{array}\right)\left(\begin{array}{c}{\mathbf{u}}_h\\ {\mathbf{v}}_h\\ p_h\end{array}\right)=\left(\begin{array}{c}\mathrm{rhs}\mathbf{u}\\ \mathrm{rhs}\mathbf{v}\\ \mathrm{rhs}p\end{array}\right),\end{array}$$

(25)

here, $\mathbf{S}$ denotes the operator for diffusion, reaction, and convection from the governing equations, whereas $\mathbf{B}$ and ${\mathbf{B}}^T$ represent the discrete gradient and divergence operators, respectively.

4. Solution Algorithm

To solve the system of nonlinear algebraic equations of saddle point type presented in (25), we employ the Newton iteration method, which can exhibit quadratic convergence for solutions that are sufficiently close. This method involves finding a root of the residual

$$\mathbf{R}\left(\mathbf{X}\right)=\mathbf{0},$$

(26)

by utilizing the already-known function value and its first derivative. The formula for the Newton iteration with damping is given as follows:

$${\mathbf{X}}^{n+1}={\mathbf{X}}^n+{\omega }^n{\left[\frac{\partial \mathbf{R}\left({\mathbf{X}}^{\mathbf{n}}\right)}{\partial \mathbf{X}}\right]}^{-1}\mathbf{R}\left({\mathbf{X}}^n\right),$$

(27)

here, $\mathbf{X}=({\mathbf{u}}_h,{\mathbf{v}}_h,p_h)$ denotes the solution vector, n represents the iteration number, and ${\omega }^n$ is the damping parameter. The Jacobian matrix $\frac{\partial \mathbf{R}\left({\mathbf{X}}^{\mathbf{n}}\right)}{\partial \mathbf{X}}$ is computed by the finite differences taking the residual matrix $\mathbf{R}\left(\mathbf{X}\right)$

$${\left[\frac{\partial \mathbf{R}\left({\mathbf{X}}^{\mathbf{n}}\right)}{\partial \mathbf{X}}\right]}_{ij}\approx \frac{{\left[\mathbf{R}\right]}_i({\mathbf{X}}^n+{\alpha }_j{\mathbf{e}}_j)-{\left[\mathbf{R}\right]}_i({\mathbf{X}}^n-{\alpha }_j{\mathbf{e}}_j)}{2{\alpha }_j},$$

(28)

In Equation (28), the coefficients ${\alpha }_j>0$ represent the increments at each iteration step n, while ${\mathbf{e}}_j$ denotes the unit basis vectors in ${\mathbf{R}}^n$. For a more detailed explanation, please refer to [26,27]. Furthermore, the parameter ${\omega }^n\in (-1,0)$ is chosen such that the error measure decreases with each iteration.

$$\mathbf{R}\left({\mathbf{X}}^{n+1}\right)·{\mathbf{X}}^{n+1}\le \mathbf{R}\left({\mathbf{X}}^n\right)·{\mathbf{X}}^n,$$

The damping parameter in the Newton iteration greatly enhances its robustness, particularly when the current approximation ${\mathbf{X}}^n$ is not close enough to the final solution. For more information, please refer to [26,28]. In this $2D$ problem, we employ a direct solver for sparse systems such as MUMPS [29]. Algorithm 1 presents the Newton iteration algorithm used in this study. Algorithm 1 outlines the Newton iteration and line search method used in this study.

Algorithm 1  Newton iteration and line search

1: Input nonlinear tolerance parameter 2: Initialize $n=0$ and take ${\mathbf{X}}^n$ as the starting guess. 3: Build the residual vector $\mathbf{R}\left({\mathbf{X}}^n\right)=\mathbf{A}{\mathbf{X}}^n-\mathbf{b}$. 4: Compute the Jacobian matrix $\mathbf{J}\left({\mathbf{X}}^n\right)=\frac{\partial \mathbf{R}}{\partial \mathbf{X}}\left({\mathbf{X}}^n\right)$. 5: Solve the linear system for correction of $\delta \mathbf{X}$: $$\mathbf{J}\left({\mathbf{X}}^n\right)\delta {\mathbf{X}}^n=\mathbf{R}\left({\mathbf{X}}^n\right)$$ 6: Find the optimal step length ${\omega }^n\in (-1,0]$. 7: Update the solution: ${\mathbf{X}}^{n+1}={\mathbf{X}}^n+{\omega }^n\delta {\mathbf{X}}^n$.

4.1. Problem Definition

The rectangular fluid channel in this problem has a length of $2.5m$ and a height of $0.4m$. The rectangular solid flap has a width of $0.04m$ and a height of $0.2m$. The complete geometry of the channel with a vertical flap is depicted in Figure 2, which includes the coarse mesh, geometrical details, and boundary conditions specification.

Table 1 gives the material properties.

For this particular problem, the channel cross-section is rectangular with a length of $2.5m$ and a height of $0.4m$. The rectangular solid flap has a width of $0.04m$ and a height of $0.2m$. A parabolic inflow velocity profile is applied on the left wall with a maximum inflow velocity of $\mathbf{v}=0.3$ given by:

$$\begin{array}{c}\hfill \mathbf{v}{\left(·,y\right)}_{inflow}=\left(\frac{1.5{\mathbf{v}}_{mean}4.y\left(0.4-y\right)}{0.4^2},0\right),\end{array}$$

here, ${\mathbf{v}}_{mean}$ is equal to $\frac{2}{3}$ of ${\mathbf{v}}_{max}$. The top and bottom walls and the vertical flap are assigned no-slip boundary conditions. On the right wall, a do-nothing boundary condition is set. Additionally, the fluid-solid interface is assumed to have equal normal stresses, i.e.,

$${\sigma }_n^f={\sigma }_n^s.$$

(29)

Figure 2 illustrates the complete geometry of the channel with a vertical flap, including the coarse mesh, geometrical details, and boundary conditions specifications. The drag coefficient and lift coefficient are computed using the following equations:

$$\begin{array}{c}\hfill C_D=\frac{2F_D}{\rho {{\mathbf{v}}_{mean}}^{2}L},C_L=\frac{2F_L}{\rho {{\mathbf{v}}_{mean}}^{2}L},\end{array}$$

where $F_D$ and $F_L$ are given by:

$$\begin{array}{c}\hfill F_D={\int }_S(\rho \eta \frac{\partial u_t}{\partial \overline{n}}{n}_y-p{n}_x)dS,F_L=-{\int }_S(\rho \eta \frac{\partial u_t}{\partial \overline{n}}{n}_x+p{n}_y)dS,\end{array}$$

respectively. Here, S represents the surface of the vertical flap, and $\overline{n}$ is the unit normal vector to S. The symbols $\rho$, p, $\eta$, and $u_t$ denote the density, pressure, dynamic viscosity, and tangential velocity, respectively. The constant L represents the characteristic length scale of the geometry, and ${\mathbf{v}}_{mean}$ is the mean inflow velocity.

5. Numerical Results

The present configuration is based on the FSI1-benchmark presented in [30], which is a well-established benchmark for validating fluid-structure interaction (FSI) simulations. The benchmark features a rectangular channel with a vertical flap, where the goal is to study the interaction between the fluid and the flexible structure. The geometry of the benchmark is relatively simple, with the channel length and height fixed at 2.5m and 0.4m, respectively. The width and height of the vertical flap are set at 0.04m and 0.2m, respectively.

The benchmark involves relatively small beam displacements (around $4\%$), making it challenging to ensure grid-independent convergence. As such, accurate numerical simulations of this benchmark require careful consideration of grid resolution, numerical scheme, and time-step size. To further validate the implementation, we consider the significant deformations of the vertical flap induced by the SVK material model, Neo-Hookean model, and Mooney-Rivlin models, which have a noticeable impact on the fluid flow behavior. Through these simulations, we aim to demonstrate the effectiveness and accuracy of our FSI solver in capturing the complex fluid-structure interaction phenomena.

5.1. Case 1: SVK material model for the vertical flap and rigid walls

In Case Section 5.1, all the outer boundaries are kept rigid while the SVK material model is applied to a vertical flap located inside the fluid channel. The figures presented in this section show the results for the SVK material model. Figure 3 provides an insight into the significant deformation of the flap due to the effect of fluid flow velocity magnitude. The surface velocity plot along with arrow representation is shown for a velocity of $\mathbf{v}=0.12m/s$. Figure 4 displays the pressure distribution within the channel and areas with high-pressure distributions. The pressure contours are presented for a velocity of $0.12m/s$. Furthermore, Figure 5 demonstrates the X and Y displacement fields of the top of the flap under the effect of fluid. The flap experiences noticeable displacement due to the applied fluid force. These results provide insight into the behavior of the SVK material model applied to the vertical flap under the effect of fluid flow velocity magnitude. The pressure distribution and displacement fields illustrate the impact of fluid flow on the flap, and these findings can be useful for further research in this field.

5.2. Case 2: SVK (Top Roof) and Neo-Hookean (Vertical Flap)

In Case Section 5.2, the SVK material model is utilized at the top solid wall, while a Neo-Hookean material model is implemented on a vertical flap within the fluid channel. Figure 6, Figure 7 and Figure 8 depict the vector magnitude, pressure distribution, and displacement on the vertical bar, providing an insight into the behavior of the system under the combined effects of fluid flow and material properties.

Figure 6 displays the surface velocity plot and arrow representation, showing the magnitude of the fluid flow and the small deformation of the flap under the fluid flow. In Figure 7, the pressure distribution and the area of great impact are clearly visible. Finally, Figure 8 illustrates the X and Y displacement fields of the top of the flap under the effect of fluid. The figures demonstrate that the Neo-Hookean material model significantly affects the deformation of the vertical flap, whereas the SVK material model has a relatively minor impact on the behavior of the top roof.

In summary, the results of Case Section 5.2 highlight the importance of selecting appropriate material models for each component of the system and the need to consider the combined effects of fluid flow and material properties in designing and analyzing fluid-structure interaction problems.

5.3. Case 3: SVK (Top Roof) and Mooney Rivlin (vertical Flap)

In Case 3, the SVK material model is assigned to the upper elastic wall, while the Mooney-Rivlin material properties are assigned to the vertical elastic flap. This case aims to understand the behavior of the Mooney-Rivlin material model under the influence of fluid flow.

The figures presented in this section, i.e., Figure 9, Figure 10 and Figure 11, provide a detailed understanding of the behavior of the Mooney-Rivlin material model under the influence of fluid flow. Figure 9 shows the vector magnitude of the fluid flow and the very little deformation of the flap under the fluid flow. The figure also shows the direction of the fluid flow using arrow representation. In Figure 10, pressure distribution, and area of great impact are quite evident. This figure highlights the maximum pressure points on the flap, providing insights into the structural integrity of the material model. Finally, Figure 11 shows the X and Y displacement fields of the top of the flap under the effect of fluid.

The results obtained from this case can be used to improve the design of materials used in various applications, such as in the aerospace industry or medical devices. By understanding the behavior of different material models under the influence of fluid flow, engineers can design more efficient and robust materials that can withstand fluid forces. Additionally, the numerical simulation approach used in this study can be extended to other complex geometries and material models to gain a better understanding of their behavior under different conditions.

5.4. Calculation of Drag and Lift

To evaluate the performance of the fluid-structure interaction, the parametric sweep function was utilized to simulate the flow at various inlet velocities. In this study, Drag and Lift forces were calculated using the localized stress components. Drag force was determined by performing fourth-order integration of stress in the x-direction, while the Lift force was computed by integrating the stress in the y-direction. These calculations were carried out by performing surface integration across the flap.

To provide a graphical representation of the Drag and Lift forces for the three cases considered, Figure 12 shows the Drag and Lift forces for Case 1, Figure 13 illustrates the Drag and Lift forces for Case 2, and Figure 14 presents the Drag and Lift forces for Case 3. These figures provide insight into the variations of Drag and Lift forces for different inlet velocities and material models.

5.5. Incorporating Non-Newtonian Fluid Model

Additionally, the Carreau fluid model has been widely used in simulating blood flow due to its accuracy in describing the non-Newtonian behavior of blood. The shear-thinning property of blood is due to the presence of large macromolecules, such as proteins and red blood cells, which can undergo deformation in response to shear stress [31]. The present study investigates the impact of non-Newtonian fluid on the hemodynamic problem at hand. To accurately model the rheological behavior of blood, the Carreau fluid model has been adopted [32]. This model accounts for the shear-thinning behavior observed in blood, as well as other non-Newtonian effects [31]. The viscosity is incorporated using the following equation:

$$\begin{array}{c}\hfill \mu \left(\right|\mathbf{D}\left(\mathbf{v}\right)\left|\right)={\mu }_{\infty }+({\mu }_0-{\mu }_{\infty })(1+K|\mathbf{D}\left(\mathbf{v}\right){{|}^2)}^n,\end{array}$$

(30)

where ${\mu }_0$ and ${\mu }_{\infty }$ are the low and high shear rate asymptotic values, $\mathbf{D}$ is the deformation gradient, and parameters n and K control the transition region. The values of the parameters are chosen to be ${\mu }_{\infty }=0.00345,{\mathrm{Ns}/\mathrm{m}}^2$, ${\mu }_0=0.056,{\mathrm{Ns}/\mathrm{m}}^2$, $K=10.976$, and $n=-0.3216$.

Figure 15 (left) shows a coarse mesh based on abstracted biomedical data, where the red and green colors represent elastic walls with separate stiffness and the blue color represents the flow region. On the right side of the same figure, a cut plane of $3D$ specimens computed tomographic image and magnetic resonance imaging data of human neurocrania is shown. The stress distribution effects at peak systole are illustrated in Figure 16. The pressure on the wall is nearly constant when the flow is fully developed, as indicated in Figure 17.

Figure 18 illustrates the behavior of blood flow toward the aneurysm, demonstrating both its magnitude and intensity. The top part of the figure displays the velocity magnitude of the entire domain, while the bottom part focuses on a zoomed-in view of the blood flow magnitude inside the aneurysm.

It is evident from the figure that the blood flow intensity is higher near the neck of the aneurysm, which could potentially lead to complications. Additionally, the flow behavior within the aneurysm is affected by the surrounding vessels, resulting in complex flow patterns.

It is worth noting that the accurate prediction of blood flow behavior is crucial for the diagnosis and treatment of aneurysms. With the advancements in computational fluid dynamics, researchers can now simulate and analyze the blood flow within aneurysms, aiding in the development of effective treatment strategies.

6. Conclusions

The coupling of solid mechanics and fluid dynamics is an important area of study in the field of biomechanics, as it allows for a better understanding of the behavior of biological tissues and organs under the influence of external forces. In this regard, the deformation of hyperelastic materials has been widely investigated due to its relevance in many biological applications, such as in the design of cardiovascular implants, tissue engineering, and soft robotics.

The findings presented above demonstrate that the deformation of a hyperelastic flap is strongly influenced by the velocity of the fluid that interacts with it. This is due to the fact that the fluid exerts forces on the surface of the flap, causing it to deform in a certain way. In addition, the study highlights the importance of the material properties of the hyperelastic flap in determining its behavior under fluid flow. Specifically, the hyperelastic properties of the material are defined by material constants such as Lame parameters and $C_{1,0}$ and $C_{0,1}$, which determine the material’s ability to resist deformation.

Understanding the behavior of hyperelastic materials under fluid flow is important for designing effective biomedical devices, such as prosthetic heart valves, stents, and arterial grafts. In these applications, the fluid flow and the deformation of the materials are intimately linked, and the design of the devices must take into account the complex interplay between the two. The findings presented above provide useful insights into this complex behavior and may help in the development of more effective biomedical devices that can withstand the complex and dynamic forces exerted by the fluid environment.

In the first case (Section 5.1), only the SVK material model was considered. The results indicate that the displacement of the flap increases with the increase in the value of the flow velocity. If the values of the Lame parameters are different, it would affect the observed displacement of the flap. This highlights the significance of material properties of the solid and fluid parameters such as velocity and media in determining the nature of the interaction between the fluid and structure. Hence, designers should carefully choose materials and adjust fluid parameters to achieve the desired outcomes for specific applications. The optimization of material properties and fluid parameters could lead to enhanced performance and efficiency in fluid-structure interaction systems.

Case (Section 5.2) introduces a more complex scenario where the fluid interacts with two hyperelastic materials. The deformation of the roof in this case is intricately linked to the deformation of the flap, resulting in a coupled response. Specifically, if the flap remains relatively undeformed as presented in the study, the flow diverts towards the roof, leading to its deformation. However, if the flap undergoes significant deformation, the diversion of flow towards the top roof is reduced, resulting in less deformation of the roof. It is worth noting that the placement of the flap also plays a significant role in determining the deformation of the roof in this scenario.

In Case (Section 5.3), two hyperelastic material models were examined in terms of their interaction with a fluid. The material properties were carefully chosen to produce minimal deformation of the flap, with most of the flow’s effects being directed towards the top roof. The displacement of the vertical flap in the y-direction was found to be significantly less than that in the x-direction. Furthermore, the velocity plot displayed the diversion of the flow. To examine the impact of velocity change on material deformation, the step size of the mean inlet velocity was decreased for Cases 2 and 3. The values of the material constants, such as $C_{1,0}$ and $C_{0,1}$, were found to play a crucial role in determining the behavior of the hyperelastic material.

The investigation presented in this study holds immense significance in the field of fluid-structure interaction, particularly concerning hyperelastic materials. The numerical approach employed in this study provides insightful information regarding the behavior of the fluid and structural variables during the interaction. Although experimentation provides accurate results, it is often challenging to conduct studies in wind tunnels to obtain deformation and displacement values of the materials due to the high cost of data acquisition and arrangement.

It is noteworthy to mention that the estimation and validation of wall shear stress (WSS) on arteries under the influence of pulsatile blood flow is of critical importance, especially in reviewing the coronary arteries for the rupture of atherosclerosis plaques. However, the lack of standard reference, inherent limitations of methodologies, and measurements call for significant improvements in accurate simulation and estimation to advance clinical research for life-saving purposes.