Search

The choice of numerical integrator in approximating solutions to dynamic partial differential equations depends on the smallest time-scale of the problem at hand. Large-scale deformations in elastic solids contain both shear waves and bulk waves, the latter of which can travel infinitely fast in incompressible materials. Explicit schemes, which are favored for their efficiency in resolving low-speed dynamics, are bound by time step size restrictions that inversely scale with the fastest wave speed. Implicit schemes can enable larger time step sizes regardless of the wave speeds present, though they are much more computationally expensive. Semi-implicit methods, which are more stable than explicit methods and more efficient than implicit methods, are emerging in the literature, though their applicability to nonlinear elasticity is not extensively studied. In this research, we develop and investigate the functionality of two time integration schemes for the resolution of large-scale dynamics in nearly- and fully-incompressible hyperelastic solids: a Modified Semi-implicit Backward Differentiation Formula integrator (MSBDF2) and a forward Euler / Semi-implicit Backward Differentiation Formula Runge-Kutta integrator (FEBDF2). We prove and empirically verify second order accuracy for both schemes. The stability properties of both methods are derived and numerically verified. We find FEBDF2 has a maximum time step size that inversely scales with the shear wave speed and is unaffected by the bulk wave speed -- the desired stability property of a semi-implicit scheme. Finally, we empirically determine that semi-implicit schemes struggle to preserve volume globally when using nonlinear incompressibility conditions, even under temporal and spatial refinement. more...

In the class of immersed boundary (IB) methods, the choice of the delta function plays a crucial role in transferring information between fluid and solid domains. Most prior work has used isotropic kernels that do not preserve the divergence-free condition of the velocity field, leading to loss of incompressibility of the solid when interpolating velocity to Lagrangian markers. To address this issue, in simulations involving large deformations of incompressible hyperelastic structures immersed in fluid, researchers often use stabilization approaches such as adding a volumetric energy term. Composite B-spline (CBS) kernels offer an alternative by maintaining the discrete divergence-free property. This work evaluates CBS kernels in terms of volume conservation and accuracy, comparing them with isotropic kernel functions using a construction introduced by Peskin (IB kernels) and B-spline (BS) kernels. Benchmark tests include pressure-loaded and shear-dominated flows, such as an elastic band under pressure loads, a pressurized membrane, a compressed block, Cook's membrane, and a slanted channel flow. Additionally, we validate our methodology using a complex fluid-structure interaction model of bioprosthetic heart valve dynamics. Results demonstrate that CBS kernels achieve superior volume conservation compared to isotropic kernels, eliminating the need for stabilization techniques. Further, CBS kernels converge on coarser fluid grids, while IB and BS kernels need finer grids for comparable accuracy. Unlike IB and BS kernels, which perform better with larger mesh ratios, CBS kernels improve with smaller mesh ratios. Wider kernels provide more accurate results across all methods, but CBS kernels are less sensitive to grid spacing variations than isotropic kernels. more...

The immersed interface method (IIM) for models of fluid flow and fluid-structure interaction imposes jump conditions that capture stress discontinuities generated by forces that are concentrated along immersed boundaries. Most prior work using the IIM for fluid dynamic applications has focused on smooth interfaces, but boundaries with sharp features such as corners and edges can appear in practical analyses, particularly on engineered structures. The present study builds on our work to integrate finite element-type representations of interface geometries with the IIM. Initial realizations of this approach used a continuous Galerkin (CG) finite element discretization for the boundary, but as we show herein, these approaches generate large errors near sharp geometrical features. To overcome this difficulty, this study introduces an IIM approach using discontinuous Galerkin (DG) representation of the jump conditions. Numerical examples explore the impacts of different interface representations on accuracy for both smooth and sharp boundaries, particularly flows interacting with fixed interface configurations. We demonstrate that using a DG approach provides accuracy that is comparable to the CG method for smooth cases. Further, we identify a time step size restriction for the CG representation that is directly related to the sharpness of the geometry. In contrast, time step size restrictions imposed by DG representations are demonstrated to be insensitive to the presence of sharp features., Comment: 21 pages, 14 figures more...

Multiphase flows are an important class of fluid flow and their study facilitates the development of diverse applications in industrial, natural and biomedical systems. Simulating such flows requires significant computational resources, making it prudent to devise an adaptive mesh refinement (AMR) method to mitigate this burden. We use a mathematical model that takes a continuum mechanical approach to describe multiphase mixture flows. The resulting system of equations poses numerical challenges due to the presence of multiple non-linear terms and a co-incompressibility condition, while the resulting fluid dynamics necessitate the development of an adaptive mesh refinement technique to accurately capture regions of interest while keeping computational costs low. We present an accurate, robust, and efficient computational method for simulating multiphase mixtures on adaptive grids, and utilize a multigrid solver to precondition the saddle-point system. We demonstrate that the AMR solver asymptotically approaches second order accuracy in $L^1$, $L^2$ and $L^\infty$ norms for all solution variables of the Newtonian and non-Newtonian models. All experiments demonstrate the solver is stable provided the time step size satisfies the imposed CFL condition. The solver can accurately resolve sharp gradients in the solution and, with the multigrid preconditioner, the solver behavior is independent of grid spacing. Our AMR solver offers a major cost savings benefit, providing up to 10x speedup in the numerical experiments presented here, with greater speedup possible depending on the problem set-up., Comment: Submitted to Journal of Computational Physics more...

This paper presents an approach to enhance volume conservation in the immersed boundary (IB) method by using regularized delta functions derived from composite B-splines. The conventional IB method, while effective for fluid-structure interaction applications, has long been challenged by poor volume conservation, particularly evident in simulations of pressurized, closed membranes. We demonstrate that composite B-spline regularized delta functions significantly enhance volume conservation through two complementary properties: they provide continuously divergence-free velocity interpolants and maintain the gradient character of forces corresponding to mean pressure jumps across interfaces. By correctly representing these forces as discrete gradients, they eliminate a key source of spurious flows that typically plague immersed boundary computations. Our approach maintains the local nature of the classical IB method, avoiding the computational overhead associated with the non-local Divergence-Free Immersed Boundary (DFIB) method's construction of an explicit velocity potential which requires additional Poisson solves for interpolation and force spreading operations. We show that sufficiently regular composite B-spline kernels maintain initial volumes to within machine precision. We provide a detailed analysis of the relationship between kernel regularity and the accuracy of force spreading and velocity interpolation operations. Our findings indicate that composite B-splines of at least $C^1$ regularity produce results comparable to the DFIB method in dynamic simulations, with errors in volume conservation dominated by truncation error of the time-stepping scheme. The proposed approach requires minimal modifications to an existing IB code, making it an accessible improvement for a wide range of applications in computational fluid dynamics and fluid-structure interaction., Comment: 35 pages, 14 figures more...

Microorganism motility often takes place within complex, viscoelastic fluid environments, e.g., sperm in cervicovaginal mucus and bacteria in biofilms. In such complex fluids, strains and stresses generated by the microorganism are stored and relax across a spectrum of length and time scales and the complex fluid can be driven out of its linear response regime. Phenomena not possible in viscous media thereby arise from feedback between the "swimmer" and the complex fluid, making swimming efficiency co-dependent on the propulsion mechanism and fluid properties. Here we parameterize a flagellar motor and filament properties together with elastic relaxation and nonlinear shear-thinning properties of the fluid in a computational immersed boundary model. We then explore swimming efficiency over this parameter space. One exemplary insight is that motor efficiency (measured by the volumetric flow rate) can be boosted vs.\ degraded by moderate vs.\ strong shear-thinning of the viscoelastic environment., Comment: 15 pages, 8 figures more...

We present and analyze a series of benchmark tests regarding the application of the immersed boundary (IB) method to viscoelastic flows through and around non-trivial, stationary geometries. The IB method is widely used for the simulation of biological fluid dynamics and other modeling scenarios where a structure is immersed in a fluid. Although the IB method has been most commonly used to model systems with viscous incompressible fluids, it also can be applied to visoelastic fluids, and has enabled the study of a wide variety of dynamical problems including the settling of vesicles and the swimming of elastic filaments in fluids modeled by the Oldroyd-B constuitive equation. However, to date, relatively little work has explored the accuracy or convergence properties of the numerical scheme. Herein, we present benchmarking results for an IB solver applied to viscoelastic flows in and around non-trivial geometries using the idealized Oldroyd-B and more realistic, polymer-entanglement-based Rolie-Poly constitutive equations. We use two-dimensional numerical test cases along with results from rheology experiments to benchmark the IB method and compare it to more complex finite element and finite volume viscoelastic flow solvers. Additionally, we analyze different choices of regularized delta function and relative Lagrangian grid spacings which allow us to identify and recommend the key choices of these numerical parameters depending on the present flow regime., Comment: 31 pages, 15 figures more...

Cardiac fluid dynamics fundamentally involves interactions between complex blood flows and the structural deformations of the muscular heart walls and the thin, flexible valve leaflets. There has been longstanding scientific, engineering, and medical interest in creating mathematical models of the heart that capture, explain, and predict these fluid-structure interactions. However, existing computational models that account for interactions among the blood, the actively contracting myocardium, and the cardiac valves are limited in their abilities to predict valve performance, resolve fine-scale flow features, or use realistic descriptions of tissue biomechanics. Here we introduce and benchmark a comprehensive mathematical model of cardiac fluid dynamics in the human heart. A unique feature of our model is that it incorporates biomechanically detailed descriptions of all major cardiac structures that are calibrated using tensile tests of human tissue specimens to reflect the heart's microstructure. Further, it is the first fluid-structure interaction model of the heart that provides anatomically and physiologically detailed representations of all four cardiac valves. We demonstrate that this integrative model generates physiologic dynamics, including realistic pressure-volume loops that automatically capture isovolumetric contraction and relaxation, and predicts fine-scale flow features. None of these outputs are prescribed; instead, they emerge from interactions within our comprehensive description of cardiac physiology. Such models can serve as tools for predicting the impacts of medical devices or clinical interventions. They also can serve as platforms for mechanistic studies of cardiac pathophysiology and dysfunction, including congenital defects, cardiomyopathies, and heart failure, that are difficult or impossible to perform in patients. more...

Fish swimming style and schooling behaviour are two important aspects of underwater locomotion. How swimming style affects fish schooling is investigated by a representative problem setup consisting of two tethered NACA0012 hydrofoils undulating in a free stream flow at various wavelengths $ \lambda^* = 0.8 - 8 $, front-back distance $ D = 0,\ 0.25,\ 0.5,\ 0.75 $, phase difference $ \phi/\pi = 0,\ 0.5,\ 1,\ 1.5 $, and lateral gap distance $ G = 0.25,\ 0.3,\ 0.35 $ with fixed Reynolds number $ Re = 5000 $, Strouhal number $ St = 0.4 $, and maximum amplitude $ A_{max} = 0.1 $. In total, 336 cases were simulated by open-source software IBAMR based on immersed boundary method. The increase in wavelength generally leads to higher thrust and lateral force, consistent with a previous study regarding a single swimmer. The highest Froude efficiency is obtained at intermediate wavelength $ \lambda^* = 2 $ for either individual or two foils as a group. At side-by-side arrangement $ D = 0 $, the thrust force upon the two foils can be equivalent, indicating a stable formation. At staggered arrangement, $ D > 0 $, the follower can take significant advantage of the leader in locomotion performance by tuning phase difference. The follower's benefit decreases with relative distance. Various combinations of wavelength and relative distance can lead to distinct flow structures, indicating a tunable stealth capacity of the schooling swimmers., Comment: 53 pages, 21 figures more...

This paper develops and benchmarks an immersed peridynamics method to simulate the deformation, damage, and failure of hyperelastic materials within a fluid-structure interaction framework. The immersed peridynamics method describes an incompressible structure immersed in a viscous incompressible fluid. It expresses the momentum equation and incompressibility constraint in Eulerian form, and it describes the structural motion and resultant forces in Lagrangian form. Coupling between Eulerian and Lagrangian variables is achieved by integral transforms with Dirac delta function kernels, as in standard immersed boundary methods. The major difference between our approach and conventional immersed boundary methods is that we use peridynamics, instead of classical continuum mechanics, to determine the structural forces. We focus on non-ordinary state-based peridynamic material descriptions that allow us to use a constitutive correspondence framework that can leverage well characterized nonlinear constitutive models of soft materials. The convergence and accuracy of our approach are compared to both conventional and immersed finite element methods using widely used benchmark problems of nonlinear incompressible elasticity. We demonstrate that the immersed peridynamics method yields comparable accuracy with similar numbers of structural degrees of freedom for several choices of the size of the peridynamic horizon. We also demonstrate that the method can generate grid-converged simulations of fluid-driven material damage growth, crack formation and propagation, and rupture under large deformations. more...

Subclinical leaflet thrombosis (SLT) is a potentially serious complication of aortic valve replacement with a bioprosthetic valve in which blood clots form on the replacement valve. SLT is associated with increased risk of transient ischemic attacks and strokes and can progress to clinical leaflet thrombosis. SLT following aortic valve replacement also may be related to subsequent structural valve deterioration, which can impair the durability of the valve replacement. Because of the difficulty in clinical imaging of SLT, models are needed to determine the mechanisms of SLT and could eventually predict which patients will develop SLT. To this end, we develop methods to simulate leaflet thrombosis that combine fluid-structure interaction and a simplified thrombosis model that allows for deposition along the moving leaflets. Additionally, this model can be adapted to model deposition or absorption along other moving boundaries. We present convergence results and quantify the model's ability to realize changes in valve opening and pressures. These new approaches are an important advancement in our tools for modeling thrombosis in which they incorporate both adhesion to the surface of the moving leaflets and feedback to the fluid-structure interaction., Comment: 29 pages, 11 figures more...

This paper introduces a sharp-interface approach to simulating fluid-structure interaction involving flexible bodies described by general nonlinear material models and across a broad range of mass density ratios. This new flexible-body immersed Lagrangian-Eulerian (ILE) approach incorporates the geometrical and domain solution flexibility of the immersed boundary (IB) method with an accuracy comparable to body-fitted approaches that sharply resolve flows and stresses up to the fluid-structure interface. Unlike many IB methods, our ILE formulation uses distinct momentum equations for the fluid and solid subregions with a Dirichlet-Neumann coupling strategy that connects fluid and solid subproblems through simple interface conditions. We use a penalty method involving two representations of the fluid-structure interface. These two representations are connected by approximate Lagrange multiplier forces that impose kinematic interface conditions. This approach also enables the use of multi-rate time stepping, which allows us to take different time step sizes for the fluid and structure subproblems. Our fluid solver relies on an immersed interface method for discrete surfaces to impose stress jump conditions along complex interfaces. The dynamics of the volumetric structural mesh are determined using a standard finite element approach to large-deformation nonlinear elasticity via a nearly incompressible solid mechanics formulation. This formulation also readily accommodates compressible structures for cases in which part of the solid boundary does not contact the incompressible fluid. Comparisons are made with computational and experimental benchmarks. We also demonstrate the capabilities of this methodology by applying it to model the transport and capture of a cylindrical blood clot in an inferior vena cava filter. more...

The immersed finite element-finite difference (IFED) method is a computational approach to modeling interactions between a fluid and an immersed structure. This method uses a finite element (FE) method to approximate the stresses and forces on a structural mesh and a finite difference (FD) method to approximate the momentum of the entire fluid-structure system on a Cartesian grid. The fundamental approach used by this method follows the immersed boundary framework for modeling fluid-structure interaction (FSI), in which a force spreading operator prolongs structural forces to a Cartesian grid, and a velocity interpolation operator restricts a velocity field defined on that grid back onto the structural mesh. Force spreading and velocity interpolation both require projecting data onto the finite element space. Consequently, evaluating either coupling operator requires solving a matrix equation at every time step. Mass lumping, in which the projection matrices are replaced by diagonal approximations, has the potential to accelerate this method considerably. Constructing the coupling operators also requires determining the locations on the structure mesh where the forces and velocities are sampled. Here we show that sampling the forces and velocities at the nodes of the structural mesh is equivalent to using lumped mass matrices in the coupling operators. A key theoretical result of our analysis is that if both of these approaches are used together, the IFED method permits the use of lumped mass matrices derived from nodal quadrature rules for any standard interpolatory element. This is different from standard FE methods, which require specialized treatments to accommodate mass lumping with higher-order shape functions. Our theoretical results are confirmed by numerical benchmarks, including standard solid mechanics tests and examination of a dynamic model of a bioprosthetic heart valve. more...

An essential ingredient of a spectral method is the choice of suitable bases for test and trial spaces. On complex domains, these bases are harder to devise, necessitating the use of domain partitioning techniques such as the spectral element method. In this study, we introduce the Projection Extension (PE) method, an approach that yields spectrally accurate solutions to various problems on complex geometries without requiring domain decomposition. This technique builds on the insights used by extension methodologies such as the immersed boundary smooth extension and Smooth Forcing Extension (SFE) methods that are designed to improve the order of accuracy of the immersed boundary method. In particular, it couples an accurate extension procedure, that functions on arbitrary domains regardless of connectedness or regularity, with a least squares minimization of the boundary conditions. The resulting procedure is stable under iterative application and straightforward to generalize to higher dimensions. Moreover, it rapidly and robustly yields exponentially convergent solutions to a number of challenging test problems, including elliptic, parabolic, Newtonian fluid flow, and viscoelastic problems., Comment: 27 pages, approved for release by Pacific Northwest National Laboratory (PNNL-SA-168413) more...

The immersed boundary (IB) method is a non-body conforming approach to fluid-structure interaction (FSI) that uses an Eulerian description of the momentum, viscosity, and incompressibility of a coupled fluid-structure system and a Lagrangian description of the deformations, stresses, and resultant forces of the immersed structure. Integral transforms with Dirac delta function kernels couple the Eulerian and Lagrangian variables, and in practice, discretizations of these integral transforms use regularized delta function kernels. Many different kernel functions have been proposed, but prior numerical work investigating the impact of the choice of kernel function on the accuracy of the methodology has been limited. This work systematically studies the effect of the choice of regularized delta function in several FSI benchmark tests using the immersed finite element/difference (IFED) method, which is an extension of the IB method that uses a finite element structural discretizations combined with a Cartesian grid finite difference method for the incompressible Navier-Stokes equations. The IFED formulation evaluates the regularized delta function on a collection of interaction points that can be chosen to be denser than the nodes of the Lagrangian mesh, and this study investigates the effect of varying the relative mesh widths of the Lagrangian and Eulerian discretizations. Our results indicate that kernels satisfying a commonly imposed even-odd condition require higher resolution to achieve similar accuracy as kernels that do not satisfy. We also find that narrower kernels are more robust and that structural meshes that are substantially coarser than the Cartesian grid can yield high accuracy for shear-dominated cases but not for cases with large normal forces. We verify our results in a large-scale FSI model of a bovine pericardial bioprosthetic heart valve in a pulse duplicator. more...

This paper presents a one-sided immersed boundary (IB) method using kernel functions constructed via a moving least squares (MLS) method. The resulting kernels effectively couple structural degrees of freedom to fluid variables on only one side of the fluid-structure interface. This reduces spurious feedback forcing and internal flows that are typically observed in IB models that use isotropic kernel functions to couple the structure to fluid degrees of freedom on both sides of the interface. The method developed here extends the original MLS methodology introduced by Vanella and Balaras (J Comput Phys, 2009). Prior IB/MLS methods have used isotropic kernel functions that coupled fluid variables on both sides of the boundary to the interfacial degrees of freedom. The original IB/MLS approach converts the cubic spline weights typically employed in MLS reconstruction into an IB kernel function that satisfies particular discrete moment conditions. This paper shows that the same approach can be used to construct one-sided kernel functions (kernel functions are referred to as generating functions in the MLS literature). We also examine the performance of the new approach for a family of kernel functions introduced by Peskin. It is demonstrated that the one-sided MLS construction tends to generate non-monotone interpolation kernels with large over- and undershoots. We present two simple weight shifting strategies to construct generating functions that are positive and monotone, which enhances the stability of the resulting IB methodology. Benchmark cases are used to test the order of accuracy and verify the one-sided IB/MLS simulations in both two and three spatial dimensions. This new IB/MLS method is also used to simulate flow over the Ahmed car model, which highlights the applicability of this methodology for modeling complex engineering flows., Comment: In Press, Journal of Computational Physics (2021) more...

We present a new discretization for advection-diffusion problems with Robin boundary conditions on complex time-dependent domains. The method is based on second order cut cell finite volume methods introduced by Bochkov et al. to discretize the Laplace operator and Robin boundary condition. To overcome the small cell problem, we use a splitting scheme that uses a semi-Lagrangian method to treat advection. We demonstrate second order accuracy in the $L^1$, $L^2$, and $L^\infty$ norms for both analytic test problems and numerical convergence studies. We also demonstrate the ability of the scheme to handle conversion of one concentration field to another across a moving boundary. more...

This paper focuses on the derivation and simulation of mathematical models describing new plasma fraction in blood for patients undergoing simultaneous extracorporeal membrane oxygenation and therapeutic plasma exchange. Models for plasma exchange with either veno-arterial or veno-venous extracorporeal membrane oxygenation are considered. Two classes of models are derived for each case, one in the form of an algebraic delay equation and another in the form of a system of delay differential equations. In special cases, our models reduce to single compartment ones for plasma exchange that have been validated with experimental data. We also show that the algebraic delay equations are forward Euler discretizations of the delay differential equations, with timesteps equal to transit times through model compartments. Numerical simulations are performed to compare different model types, to investigate the impact of plasma device port switching on the efficiency of the exchange process, and to study the sensitivity of the models to their parameters. more...

An immersed interface-lattice Boltzmann method (II-LBM) is developed for modelling fluid-structure systems. The key element of this approach is the determination of the jump conditions that are satisfied by the distribution functions within the framework of the lattice Boltzmann method when forces are imposed along a surface immersed in an incompressible fluid. In this initial II-LBM, the discontinuity related to the normal portion of the interfacial force is sharply resolved by imposing the relevant jump conditions using an approach that is analogous to imposing the corresponding pressure jump condition in the incompressible Navier-Stokes equations. We show that the jump conditions for the distribution functions are the same in both single-relaxation-time and multi-relaxation-time LBM formulations. Tangential forces are treated using the immersed boundary-lattice Boltzmann method (IB-LBM). The performance of the II-LBM method is compared to both the direct forcing IB-LBM for rigid-body fluid-structure interaction, and the classical IB-LBM for elastic interfaces. Higher order accuracy is observed with the II-LBM as compared to the IB-LBM for selected benchmark problems. Because the jump conditions of the distribution function also satisfy the continuity of the velocity field across the interface, the error in the velocity field is much smaller for the II-LBM than the IB-LBM. The II-LBM is also demonstrated to provide superior volume conservation when simulating flexible boundaries. more...

This paper introduces a sharp interface method to simulate fluid-structure interaction (FSI) involving rigid bodies immersed in viscous incompressible fluids. The capabilities of this methodology are demonstrated for a range of benchmark test cases along with large-scale models of biomedical FSI. The numerical approach developed herein, which we refer to as an immersed Lagrangian-Eulerian method, integrates aspects of partitioned and immersed FSI formulations by solving separate momentum equations for the fluid and solid subdomains, as in a partitioned formulation, while also using non-conforming discretizations of the dynamic fluid and structure regions, as in an immersed formulation. A Dirichlet-Neumann coupling scheme is used, in which the motion of the immersed solid is driven by fluid traction forces evaluated along the fluid-structure interface, and the motion of the fluid along that interface is constrained to match the solid velocity and thereby satisfy the no-slip condition. To develop a practical numerical method, we adopt a penalty approach that approximately imposes the no-slip condition along the fluid-structure interface. Our fluid-structure interaction scheme relies on an immersed interface method for discrete geometries, which enables the accurate determination of both velocities and stresses along complex fluid-structure interfaces. Unlike commonly used partitioned FSI methods, which can suffer from so-called added mass effect instabilities, our methodology retains stability for test cases involving extremely small, nearly equal, equal, and large solid-fluid density ratios without requiring subiterations or complex handling of the pressure. Biomedical FSI demonstration cases are also presented including the dynamics of a bileaflet mechanical heart valve in a pulse duplicator, and transport of blood clots in a patient-averaged anatomical model of the inferior vena cava. more...

High-order numerical methods for solving elliptic equations over arbitrary domains typically require specialized machinery, such as high-quality conforming grids for finite elements method, and quadrature rules for boundary integral methods. These tools make it difficult to apply these techniques to higher dimensions. In contrast, fixed Cartesian grid methods, such as the immersed boundary (IB) method, are easy to apply and generalize, but typically are low-order accurate. In this study, we introduce the Smooth Forcing Extension (SFE) method, a fixed Cartesian grid technique that builds on the insights of the IB method, and allows one to obtain arbitrary orders of accuracy. Our approach relies on a novel Fourier continuation method to compute extensions of the inhomogeneous terms to any desired regularity. This is combined with the highly accurate Non-Uniform Fast Fourier Transform for interpolation operations to yield a fast and robust method. Numerical tests confirm that the technique performs precisely as expected on one-dimensional test problems. In higher dimensions, the performance is even better, in some cases yielding sub-geometric convergence. We also demonstrate how this technique can be applied to solving parabolic problems and for computing the eigenvalues of elliptic operators on general domains, in the process illustrating its stability and amenability to generalization. more...

The immersed boundary-finite element method (IBFE) is an approach to describing the dynamics of an elastic structure immersed in an incompressible viscous fluid. In this formulation, there are discontinuities in the pressure and viscous stress at fluid-structure interfaces. The standard immersed boundary approach, which connects the Lagrangian and Eulerian variables via integral transforms with regularized Dirac delta function kernels, smooths out these discontinuities, which generally leads to low order accuracy. This paper describes an approach to accurately resolve pressure discontinuities for these types of formulations, in which the solid may undergo large deformations. Our strategy is to decompose the physical pressure field into a sum of two pressure-like fields, one defined on the entire computational domain, which includes both the fluid and solid subregions, and one defined only on the solid subregion. Each of these fields is continuous on its domain of definition, which enables high accuracy via standard discretization methods without sacrificing sharp resolution of the pressure discontinuity. Numerical tests demonstrate that this method improves rates of convergence for displacements, velocities, stresses, and pressures, as compared to the conventional IBFE method. Further, it produces much smaller errors at reasonable numbers of degrees of freedom. The performance of this method is tested on several cases with analytic solutions, a nontrivial benchmark problem of incompressible solid mechanics, and an example involving a thick, actively contracting torus. more...

Fluid-structure systems occur in a range of scientific and engineering applications. The immersed boundary(IB) method is a widely recognized and effective modeling paradigm for simulating fluid-structure interaction(FSI) in such systems, but a difficulty of the IB formulation is that the pressure and viscous stress are generally discontinuous at the interface. The conventional IB method regularizes these discontinuities, which typically yields low-order accuracy at these interfaces. The immersed interface method(IIM) is an IB-like approach to FSI that sharply imposes stress jump conditions, enabling higher-order accuracy, but prior applications of the IIM have been largely restricted to methods that rely on smooth representations of the interface geometry. This paper introduces an IIM that uses only a C0 representation of the interface,such as those provided by standard nodal Lagrangian FE methods. Verification examples for models with prescribed motion demonstrate that the method sharply resolves stress discontinuities along the IB while avoiding the need for analytic information of the interface geometry. We demonstrate that only the lowest-order jump conditions for the pressure and velocity gradient are required to realize global 2nd-order accuracy. Specifically,we show 2nd-order global convergence rate along with nearly 2nd-order local convergence in the Eulerian velocity, and between 1st-and 2nd-order global convergence rates along with 1st-order local convergence for the Eulerian pressure. We also show 2nd-order local convergence in the interfacial displacement and velocity along with 1st-order local convergence in the fluid traction. As a demonstration of the method's ability to tackle complex geometries,this approach is also used to simulate flow in an anatomical model of the inferior vena cava., Comment: - Added a non-axisymmetric example (flow within eccentric rotating cylinder in Sec. 4.3) - Added a more in-depth analysis and comparison with a body-fitted approach for the application in Sec. 4.7 more...

The immersed boundary method is a mathematical framework for modeling fluid-structure interaction. This formulation describes the momentum, viscosity, and incompressibility of the fluid-structure system in Eulerian form, and it uses Lagrangian coordinates to describe the structural deformations, stresses, and resultant forces. Integral transforms with Dirac delta function kernels connect the Eulerian and Lagrangian frames. The fluid and the structure are both typically treated as incompressible materials. Upon discretization, however, the incompressibility of the structure is only maintained approximately. To obtain an immersed method for incompressible hyperelastic structures that is robust under large structural deformations, we introduce a volumetric energy in the solid region that stabilizes the formulation and improves the accuracy of the numerical scheme. This formulation augments the discrete Lagrange multiplier for the incompressibility constraint, thereby improving the original method's accuracy. This volumetric energy is incorporated by decomposing the strain energy into isochoric and dilatational components, as in standard solid mechanics formulations of nearly incompressible elasticity. We study the performance of the stabilized method using several quasi-static solid mechanics benchmarks, a dynamic fluid-structure interaction benchmark, and a detailed three-dimensional model of esophageal transport. The accuracy achieved by the stabilized immersed formulation is comparable to that of a stabilized finite element method for incompressible elasticity using similar numbers of structural degrees of freedom., Comment: Extensive revisions more...

This paper presents a robust, adaptive numerical scheme for simulating high density ratio and high shear multiphase flows on locally refined Cartesian grids that adapt to the evolving interfaces and track regions of high vorticity. The algorithm combines the interface capturing level set method with a variable-coefficient incompressible Navier-Stokes solver that is demonstrated to stably resolve material contrast ratios of up to six orders of magnitude. The discretization approach ensures second-order pointwise accuracy for both velocity and pressure with several physical boundary treatments, including velocity and traction boundary conditions. The paper includes several test cases that demonstrate the order of accuracy and algorithmic scalability of the flow solver. To ensure the stability of the numerical scheme in the presence of high density and viscosity ratios, we employ a consistent treatment of mass and momentum transport in the conservative form of discrete equations. This consistency is achieved by solving an additional mass balance equation, which we approximate via a strong stability preserving Runga-Kutta time integrator and by employing the same mass flux (obtained from the mass equation) in the discrete momentum equation. The scheme uses higher-order total variation diminishing (TVD) and convection-boundedness criterion (CBC) satisfying limiter to avoid numerical fluctuations in the transported density field. The high-order bounded convective transport is done on a dimension-by-dimension basis, which makes the scheme simple to implement. We also demonstrate through several test cases that the lack of consistent mass and momentum transport in non-conservative formulations, which are commonly used in practice, or the use of non-CBC satisfying limiters can yield very large numerical error and very poor accuracy for convection-dominant high density ratio flows., Comment: Figures are compressed to comply with arXiv size requirements more...

We propose a reciprocal, self-propelled model swimmer at intermediate Reynolds numbers ($Re$). Our swimmer consists of two unequal spheres that oscillate in antiphase generating nonlinear steady streaming (SS) flows. We show computationally that the SS flows enable the swimmer to propel itself, and also switch direction as $Re$ increases. We quantify the transition in the swimming direction by collapsing our data on a critical $Re$ and show that the transition in swimming directions corresponds to the reversal of the SS flows. Based on our findings, we propose that SS can be an important physical mechanism for motility at intermediate $Re$., Comment: 6 pages, 3 figures more...

In standard models of cardiac electrophysiology, including the bidomain and monodomain models, local perturbations can propagate at infinite speed. We address this unrealistic property by developing a hyperbolic bidomain model that is based on a generalization of Ohm's law with a Cattaneo-type model for the fluxes. Further, we obtain a hyperbolic monodomain model in the case that the intracellular and extracellular conductivity tensors have the same anisotropy ratio. In one spatial dimension, the hyperbolic monodomain model is equivalent to a cable model that includes axial inductances, and the relaxation times of the Cattaneo fluxes are strictly related to these inductances. A purely linear analysis shows that the inductances are negligible, but models of cardiac electrophysiology are highly nonlinear, and linear predictions may not capture the fully nonlinear dynamics. In fact, contrary to the linear analysis, we show that for simple nonlinear ionic models, an increase in conduction velocity is obtained for small and moderate values of the relaxation time. A similar behavior is also demonstrated with biophysically detailed ionic models. Using the Fenton-Karma model along with a low-order finite element spatial discretization, we numerically analyze differences between the standard monodomain model and the hyperbolic monodomain model. In a simple benchmark test, we show that the propagation of the action potential is strongly influenced by the alignment of the fibers with respect to the mesh in both the parabolic and hyperbolic models when using relatively coarse spatial discretizations. Accurate predictions of the conduction velocity require computational mesh spacings on the order of a single cardiac cell. We also compare the two formulations in the case of spiral break up and atrial fibrillation in an anatomically detailed model of the left atrium, and [...]., Comment: 20 pages, 12 figures more...

Each year, approximately 300,000 heart valve repair or replacement procedures are performed worldwide, including approximately 70,000 aortic valve replacement surgeries in the United States alone. This paper describes progress in constructing anatomically and physiologically realistic immersed boundary (IB) models of the dynamics of the aortic root and ascending aorta. This work builds on earlier IB models of fluid-structure interaction (FSI) in the aortic root, which previously achieved realistic hemodynamics over multiple cardiac cycles, but which also were limited to simplified aortic geometries and idealized descriptions of the biomechanics of the aortic valve cusps. By contrast, the model described herein uses an anatomical geometry reconstructed from patient-specific computed tomography angiography (CTA) data, and employs a description of the elasticity of the aortic valve leaflets based on a fiber-reinforced constitutive model fit to experimental tensile test data. Numerical tests show that the model is able to resolve the leaflet biomechanics in diastole and early systole at practical grid spacings. The model is also used to examine differences in the mechanics and fluid dynamics yielded by fresh valve leaflets and glutaraldehyde-fixed leaflets similar to those used in bioprosthetic heart valves. Although there are large differences in the leaflet deformations during diastole, the differences in the open configurations of the valve models are relatively small, and nearly identical hemodynamics are obtained in all cases considered. more...

The immersed boundary (IB) method is a mathematical and numerical framework for problems of fluid-structure interaction, treating the particular case in which an elastic structure is immersed in a viscous incompressible fluid. The IB approach to such problems is to describe the elasticity of the immersed structure in Lagrangian form, and to describe the momentum, viscosity, and incompressibility of the coupled fluid-structure system in Eulerian form. Interaction between Lagrangian and Eulerian variables is mediated by integral equations with Dirac delta function kernels. The IB method provides a unified formulation for fluid-structure interaction models involving both thin elastic boundaries and also thick viscoelastic bodies. In this work, we describe the application of an adaptive, staggered-grid version of the IB method to the three-dimensional simulation of the fluid dynamics of the aortic heart valve. Our model describes the thin leaflets of the aortic valve as immersed elastic boundaries, and describes the wall of the aortic root as a thick, semi-rigid elastic structure. A physiological left-ventricular pressure waveform is used to drive flow through the model valve, and dynamic pressure loading conditions are provided by a reduced (zero-dimensional) circulation model that has been fit to clinical data. We use this model and method to simulate aortic valve dynamics over multiple cardiac cycles. The model is shown to approach rapidly a periodic steady state in which physiological cardiac output is obtained at physiological pressures. These realistic flow rates are not specified in the model, however. Instead, they emerge from the fluid-structure interaction simulation. more...

The Immersed Boundary (IB) method is a mathematical framework for constructing robust numerical methods to study fluid-structure interaction in problems involving an elastic structure immersed in a viscous fluid. The IB formulation uses an Eulerian representation of the fluid and a Lagrangian representation of the structure. The Lagrangian and Eulerian frames are coupled by integral transforms with delta function kernels. The discretized IB equations use approximations to these transforms with regularized delta function kernels to interpolate the fluid velocity to the structure, and to spread structural forces to the fluid. It is well-known that the conventional IB method can suffer from poor volume conservation since the interpolated Lagrangian velocity field is not generally divergence-free, and so this can cause spurious volume changes. In practice, the lack of volume conservation is especially pronounced for cases where there are large pressure differences across thin structural boundaries. The aim of this paper is to greatly reduce the volume error of the IB method by introducing velocity-interpolation and force-spreading schemes with the properties that the interpolated velocity field in which the structure moves is at least C1 and satisfies a continuous divergence-free condition, and that the force-spreading operator is the adjoint of the velocity-interpolation operator. We confirm through numerical experiments in two and three spatial dimensions that this new IB method is able to achieve substantial improvement in volume conservation compared to other existing IB methods, at the expense of a modest increase in the computational cost. Further, the new method provides smoother Lagrangian forces (tractions) than traditional IB methods. The method presented here is restricted to periodic computational domains. Its generalization to non-periodic domains is important future work., Comment: 49 pages, 13 figures more...

The immersed boundary method is an approach to fluid-structure interaction that uses a Lagrangian description of the structural deformations, stresses, and forces along with an Eulerian description of the momentum, viscosity, and incompressibility of the fluid-structure system. The original immersed boundary methods described immersed elastic structures using systems of flexible fibers, and even now, most immersed boundary methods still require Lagrangian meshes that are finer than the Eulerian grid. This work introduces a coupling scheme for the immersed boundary method to link the Lagrangian and Eulerian variables that facilitates independent spatial discretizations for the structure and background grid. This approach employs a finite element discretization of the structure while retaining a finite difference scheme for the Eulerian variables. We apply this method to benchmark problems involving elastic, rigid, and actively contracting structures, including an idealized model of the left ventricle of the heart. Our tests include cases in which, for a fixed Eulerian grid spacing, coarser Lagrangian structural meshes yield discretization errors that are as much as several orders of magnitude smaller than errors obtained using finer structural meshes. The Lagrangian-Eulerian coupling approach developed in this work enables the effective use of these coarse structural meshes with the immersed boundary method. This work also contrasts two different weak forms of the equations, one of which is demonstrated to be more effective for the coarse structural discretizations facilitated by our coupling approach. more...

The immersed boundary (IB) method is a widely used approach to simulating fluid-structure interaction (FSI). Although explicit versions of the IB method can suffer from severe time step size restrictions, these methods remain popular because of their simplicity and generality. In prior work (Guy et al., Adv Comput Math, 2015), some of us developed a geometric multigrid preconditioner for a stable semi-implicit IB method under Stokes flow conditions; however, this solver methodology used a Vanka-type smoother that presented limited opportunities for parallelization. This work extends this Stokes-IB solver methodology by developing smoothing techniques that are suitable for parallel implementation. Specifically, we demonstrate that an additive version of the Vanka smoother can yield an effective multigrid preconditioner for the Stokes-IB equations, and we introduce an efficient Schur complement-based smoother that is also shown to be effective for the Stokes-IB equations. We investigate the performance of these solvers for a broad range of material stiffnesses, both for Stokes flows and flows at nonzero Reynolds numbers, and for thick and thin structural models. We show here that linear solver performance degrades with increasing Reynolds number and material stiffness, especially for thin interface cases. Nonetheless, the proposed approaches promise to yield effective solution algorithms, especially at lower Reynolds numbers and at modest-to-high elastic stiffnesses. more...

In this work, we extend our previous esophageal transport model using an immersed boundary (IB) method with discrete fiber-based structures, to one using a continuum mechanics-based model that is approximated based on finite elements (IB-FE). To deal with the leakage of flow when the Lagrangian mesh becomes coarser than the fluid mesh, we employ adaptive interaction quadrature points for Lagrangian-Eulerian interaction equations based on a previous work. In particular, we introduce a new anisotropic adaptive interaction quadrature rule. The new rule permits us to vary the interaction quadrature points not only at each time-step and element but also at different orientations per element. For the material model, we extend our previous fiber-based model to a continuum-based model. We first study a case in which a three-dimensional short tube is dilated. Results match very well with those from the implicit FE method. We remark that in our IB-FE case, the three-dimensional tube undergoes a very large deformation and the Lagrangian mesh-size becomes about 6 times of Eulerian mesh-size. To validate the method in handling fiber-matrix material models, we perform a second study on dilating a long fiber-reinforced tube. Errors are small when we compare numerical solutions with analytical solutions. The technique is then applied to the problem of esophageal transport. We present three cases that differ in the material model and muscle fiber architecture. The overall transport features are consistent with those from the previous model. We remark that the continuum-based model can handle more realistic and complicated material behavior. This is demonstrated in our third case with spatially varying fiber architecture. We find this unique muscle fiber architecture could generate a so-called pressure transition zone. This suggests an important role of muscle fiber architecture in esophageal transport., Comment: 47 pages, 17 figures more...