Freeman Scholar Lecture (2021)—Sharp-Interface Immersed Boundary Methods in Fluid Dynamics

Fifty years have elapsed since Charles Peskin introduced the immersed boundary method (IBM) in his Ph.D. thesis [1], a method that allowed flow simulations around complex moving boundaries to be implemented on body-nonconformal Cartesian grids. From its origin as a niche tool for biologists and applied mathematicians, IBMs have now become widely applied in various domains of fluid dynamics. Early development, up to the 1990s, centered on “diffuse-interface” methods [2–7]. However, the last 25 years have witnessed significant advancements in “sharp-interface” IBMs, which form the focus of this article. We will present here an overview of sharp-interface IB methods developed in our group, accompanied by a diverse set of applications showcasing their versatility.

Before delving into the specifics of the sharp-interface methods addressed in this paper, it is important to establish definitions for “diffuse” and “sharp” interface IB methods. These terms find their roots in the physics of fluid interfaces [7,8], and their use in categorizing IB methods is not as obvious as it seems. In our view, when an IB method is categorized as a “sharp-interface method,” it implies that the method imposes boundary conditions on the immersed boundary with a level of accuracy and precision comparable to an equivalent body-conformable grid (BCG) method. This equivalence to BCG methods (see Fig. 1) can be summarized by the following set of conditions for sharp-interface IBMs [9], which have been fundamental in guiding our IBM development efforts:

1. The imposition of no-slip, no-penetration boundary conditions is at a set of discrete points located precisely on the immersed boundary and nowhere else within the flow domain.

2. The spacing and resolution of these boundary points align with the underlying fluid grid.

3. The boundary conditions on these boundary points are imposed with an accuracy that is consistent with the underlying numerical scheme employed.

4. All grid points/cells in the flow domain impose the “native” governing equation, i.e., the Navier–Stokes equation, and not any other augmented or auxiliary equation.

5. The governing equations of flow are not solved inside the immersed body.

A necessary condition for an IB method to be considered a sharp-interface method is that for an immersed boundary that aligns perfectly with the underlying grid, the IB method should be entirely consistent with the corresponding body-conformal grid method. However, this condition is necessary but not sufficient. For instance, the “stair-step” or “staircase” method [10] represents the immersed boundary with a stair-step shaped surface that aligns with the underlying Cartesian grid. This method satisfies the above BCG consistency condition but does not meet the first condition in the itemized list above and therefore does not classify as a sharp-interface method.

Peskin's method [1], along with feedback forcing [4], fictitious domain [11], and front tracking methods [12], incorporates a localized body-force term in the continuous equations to impose the boundary condition on the body surface. However, for implementation in the discrete form of the equations, this localized force in the continuous equations is regularized, i.e., spread across a finite number of grid points. These methods therefore fall under the category of diffuse-interface methods because they do not satisfy most of the five conditions. Other approaches fail to satisfy one or more of these five criteria. For example, the mask method [13] distributes the discrete force across a layer of grid cells, qualifying it as a diffuse-interface technique.

We also emphasize that the conditions for sharpness enumerated above may neither be necessary nor sufficient for generating high-quality results from IBM simulations. This is because factors such as the choice of the spatial and temporal discretization schemes [14] may be entangled with any particular implementation of IBM in a way as to impose additional constraints on its accuracy and fidelity. For instance, some newer methods such as those in Refs. [15] and [16] do not satisfy one or more of the conditions enumerated above, but nevertheless provide good accuracy and fidelity. Thus, the above set of conditions are best treated as a general rubric for categorizing any given IBM and not necessarily for a direct assessment of the quality of the simulations. In Secs. 2.1 and 2.2, we describe two sharp-interface IB methods that have been developed and employed by our group.

In cut-cell methods, commonly used with finite volume discretizations, the cells intersected by the immersed boundary are “reshaped” to fit the local boundary (see Fig. 2), and their discretization is accordingly adjusted. Applying cut-cell methods to three-dimensional (3D) problems is challenging and infrequent due to the complexity of dealing with cut-cells of various topologies [18,19]. To address this challenge, Seo and Mittal [20] proposed a “virtual” cut-cell method, which, when combined with a 3D finite difference ghost-cell IBM [21], achieved a higher level of discrete conservation while avoiding the difficulties associated with 3D cut-cells.

Cut-cell methods [17,19,22,23] are unequivocally sharp-interface methods since they clearly satisfy the conditions listed earlier. In fact, the cut-cell method can be viewed as a body-conformal finite volume method where the finite volumes are Cartesian everywhere in the domain except at the boundary, where they assume non-Cartesian shapes (see Fig. 1(b)).

Within the category of discrete forcing methods [24] are methods that employ “ghost-cells” (see Fig. 3) to impose the boundary conditions on the immersed boundary. Using a layer of cells just outside the computational boundary to apply external boundary conditions is a common practice in many BCG codes [25–27]. This approach has also been utilized to apply boundary conditions on immersed no-slip boundaries in body-nonconformal grid codes [20,21,28–30] and for multifluid interfaces [31,32]. Typically employed in finite difference methods, these ghost-cell techniques do not reshape boundary cells. Instead, they use an auxiliary equation for the ghost-cells to enforce the boundary condition on the immersed boundary, which is coupled to the flow variables near the boundary. In our own solver [21], the value of the ghost-cell (GC) velocity is determined through a second-order centered approximation, using the boundary condition at the body-intercept (BI) point and the flow variable interpolated at the image-point (IP). The pressure gradient boundary condition, which employs the same approximation, is formally first-order accurate, but this results in a globally second-order accurate approximation for pressure.

The ghost-cell method can also be extended to higher orders. In the work of Seo and Mittal [33], the ghost-cell value is obtained from an nth order approximating polynomial interpolation by using an arbitrarily large number of stencil points around the immersed boundary (see Fig. 4). This high-order ghost-cell method has been employed for the wave propagation [33,34] and compressible flow problems [35,36].

The ghost-cell methods by Tseng and Ferziger [28], Ghias et al. [30], Mittal et al. [21], and Seo and Mittal [20], as well as embedded boundary methods [37–40], immersed interface methods [18], and the level-set based Cartesian grid methods [41,42], meet the conditions for sharp IB methods. Many of these methods enforce velocity boundary conditions on the immersed boundary with local second-order accuracy which exceeds the requirement for global second-order accuracy, thereby offering high resolution to the boundary layers on the immersed surfaces. The sharp-interface method by Seo and Mittal [20] even integrates discrete conservation of mass in localized regions near the immersed boundary through a virtual cut-cell method, enhancing the accuracy of boundary layer flows. Finally, unlike cut-cell methods, ghost-cell methods are more easily extended to 3D problems, and many of the above references showcase applications of these methods to a variety of 3D problems with complex immersed boundaries. Consequently, the ghost-cell based IBM described above [20,21] and the resulting flow solver (called “vicar3d”) has become our go-to solver, and all the incompressible flow simulations described later in the paper have been carried out with this solver.

The fifth condition in the list above may seem harmless initially, but it has significant implications for IB methods. This condition relates to the lack of an explicit pressure boundary condition on the immersed surface in certain IB methods, such as continuous forcing and most penalization techniques. Consequently, these methods require the computation of pressure and velocity throughout the entire grid, including within the immersed body. This approach can greatly reduce the computational cost of the simulation, since imposing a pressure boundary condition on the immersed surface could increase the computational effort for solving the pressure Poisson equation (PPE) by an order of magnitude or more.

As noted in Ref. [9], in the context of the fractional-step method, commonly used by these solvers, the slip and penetration on a surface at the end of the pressure correction step are $O[(Δt)m∂p/∂τ]$ and $O[(Δt)m∂p/∂n]$⁠, respectively, where m is the order of the fractional-step method and $τ$ and n are the directions tangential and normal to the immersed body. For standard fractional-step schemes, $m=1$⁠, but second-order implementations [43] also exist. In BCG methods, it is standard to apply $∂p/∂n=0$ on all boundaries since it is consistent with the fractional-step method [44,45] and this ensures satisfaction of no-penetration to machine zero at the end of the time-step irrespective of the size of the time-step $Δt$ as well as the magnitude of the normal pressure gradient. The latter can be large, especially in regions where the flow impacts the body normally. The cut-cell and ghost-cell methods developed by our group [17,21] and others [41,42] explicitly impose the Neumann condition for pressure on the immersed boundary. This allows us to enforce no-penetration to machine zero at the boundary points at the end of the time-step, decoupling the grid points inside the immersed body from those outside and eliminating the need to solve for the flow variables inside the body.

Cut-cell methods, which rely on finite volume discretization, ensure discrete conservation of mass and momentum. However, as far as we are aware, all other IB methods inherently fail to achieve discrete conservation for cells near the immersed boundary since discrete conservation is not embedded in the immersed boundary formulations of any IB method. As with body-conformal methods, lack of discrete conservation is particularly important when the grid resolution is marginal. In sharp IB methods, this lack of discrete conservation can lead to undesirable effects like spurious temporal oscillations in surface quantities [20]. To address geometric volume conservation and enhance mass conservation, we have implemented a virtual cut-cell method in our sharp-interface IBM solver [20]. This approach uses cut-cell volumes to discretize the PPE, improving volume and mass conservation. The momentum equation is solved using a second-order central finite difference discretization with the ghost-cell method [21], eliminating the need for additional computations of momentum fluxes and stresses for the cut-cells. Cells intersected by the immersed boundary are modified as shown in Fig. 5. For regular cut-cells, the PPE is discretized using a finite volume method. The small cells are treated by the virtual cell-merging technique [20], in which merging is effected by transferring the source term of the Poisson equation for the small-cells to the adjacent regular and noncut cells. The value of the right-hand side of the PPE is computed for all cells including the small-cells, and then, for the small cells, the computed value is transferred to adjoining nonsmall cells (regular cut cells or normal, noncut cells) which share a face with the small-cut cell as shown in Fig. 5. This method has been shown to retain all the desirable properties of the original sharp-interface IBM, while at the same time, reducing pressure oscillations for moving boundaries by improving the geometrical volume and mass conservation. This treatment does slow down the iterative convergence of the pressure Poisson equation, but in our view, this cost is well worth the advantages in accuracy and robustness bestowed on the computer solution. For membraneous bodies (such as the bat wing and fish fin showed later), the virtual cut-cell merging can be implemented everywhere on the immersed surface except at the edges of the membrane. The virtual cut-cell can also be implemented for cases with multiple bodies (such as the fish school shown later) as long as there is a gap of a few grid cells between the two body surfaces.

For the ghost-cell method with higher accuracy, a technique that uses high-order polynomial interpolation combined with weighted least-squares error minimization has been proposed [33]. In this method, the ghost point value is determined by enforcing the boundary condition at the body-intercept point using a high-order polynomial. Specifically, a general variable $ϕ$ is approximated near the body-intercept point $(BI)$ (Fig. 4) with an Nth-degree polynomial. Fluid data points surrounding the body-intercept point are gathered within a circular (or spherical in 3D) region of radius R to solve for the unknown coefficients, which are then obtained by minimizing the least-square error. This approach has demonstrated a significant reduction in numerical dissipation and dispersion errors near the immersed boundary when using third or fourth-order polynomials [33]. The method has been successfully applied in areas such as aeroacoustics [33], wave propagation in biological tissues [34], and transonic as well as supersonic flows [36].

In this section, we summarize a range of applications of the IB methods developed in our group. These applications range from biolocomotion and physiology to compressible aerodynamics and aeroacoustic noise.

Bats are the only mammals capable of powered flight, and they employ an elastic membrane “hand-wing” for flight. The hand-wing comprises of an anisotropic membrane supported by a skeletal structure which is similar to a human hand but with highly elongated bones [46]. The hand-wing is actuated through muscles, and the membrane then interacts with the surrounding air to generate the aerodynamic forces of lift and thrust to power flight. The other distinct feature is the high degree of articulation in the hand-wing, which provides the bat with precision control over its flapping wing kinematics, thereby allowing it to perform complex maneuvers. The interaction of airflow with the soft-membrane wing presents a complex fluid–structure interaction (FSI) problem, an understanding of which could reveal the secrets of the flight capabilities of bats and help engineers design novel bio-inspired engineered wings and aerial vehicles [47].

Here, we present our effort to simulate the complex FSI physics of a bat's membrane wing in flight [48,49]. This requires a coupling between computational fluid dynamics and computational structural dynamics, and due to large wing strokes and high deformation in the wing, IB methods are the well-suited for performing these simulations. For the flow simulations, we employ the sharp-interface IB solver vicar3d [21]. For the dynamics of the elastic wing, we have developed a new computational structural dynamics solver based on a spring network model [39,50], and the arrangement of the springs is shown in Fig. 6. We utilized the joint kinematics of the skeleton from experiments [51] and used them to prescribe the motion of the skeleton in our simulations. The material properties of the wing were obtained from Li et al. [52]. The actual Reynolds numbers (defined here as $U∞A/ν$⁠) for bat flight could range from $104$ to $105$⁠, but given the expense of resolving such flows, we choose a Reynolds number of 1500, which is high enough to generate a complex vortex-dominated flow that would be representative of the flow at higher Reynolds numbers, but low enough to enable a resolved simulation. The Strouhal number (⁠$fA/U∞$⁠) based on freestream velocity (⁠$U∞$⁠) and stroke amplitude (A) for the simulation was chosen to be 1.0. This value is on the higher end within the flight regime observed for bats [53] and is expected to elicit complex flow and structural dynamics in the membrane, which would challenge the computational method.

Figure 7 shows results of the coupled FSI simulations for a model based on Cynopterus brachyotis. The snapshots at $t/T=0−0.3$ correspond to the upstroke, and snapshots at $t/T=0.45−0.75$ simulation results during the downstroke. Figure 8 shows the time-varying lift and thrust forces generated by the flapping motion of the wing. The aerodynamic forces show oscillatory behavior during the flapping cycle, and this could either be due to high-frequency vortex shedding or a flutter in the wing membrane. Due to the coupled nature of this simulation, dissecting the contribution from these two effects is difficult, and here we employ the force partitioning method [54] to decompose the lift and thrust forces into their contributions from vortex-induced and added-mass-induced effects, as shown in Fig. 8. It can be seen that the oscillations are predominantly due to added-mass effects associated with the flutter in the membrane.

In our past work [55], we have examined carangiform swimmers (in isolation and in pairs) at a Reynolds number (⁠$ReL=L2f/ν$⁠, where L is the body length and f is the tail-beat frequency) of 5000. Based on the size of fish, however, the Reynolds number can be much higher, and it is of interest to understand how the wake dynamics change with this important parameter. The ghost-cell based sharp-interface immersed boundary method used here [21] is capable of resolving complex flow structures generated by the swimming fish at these relatively high Reynolds numbers. These simulations were carried out for a “tethered” fish in an incoming uniform flow, and for each choice of Reynolds number, a series of simulations were carried out to determine an incoming flow velocity that results in a nearly net zero hydrodynamic force in the surge.

Figures 9(a) and 9(b) show the vortical structures of swimming fish at $ReL=5000$ and 50,000. The simulations are carried out on high-resolution computational grids employing about 233 × 10⁶ mesh points, and it was confirmed that the key features of the flow and the hydrodynamic forces were insensitive to the grid, particularly for the higher Reynolds number simulation. One can see very complex vortical structures in the wake of $ReL=50,000$ case as the vortices break down to smaller structures at this high Reynolds number. The instantaneous velocity profiles showing the boundary layer around the fish body are plotted in Fig. 9(c) for both $ReL=5000$ and 50,000 cases. The thinner boundary layer for the higher Reynolds number case is shown to be well resolved by the sharp-interface IBM. At $ReL=5000$⁠, the Strouhal number based on the tail-beat amplitude resulting in free-swimming is found to be $StA=fA/U=0.42$⁠, where A is the peak-to-peak amplitude and U is the swimming speed. At the higher Reynolds number, $ReL=50,000$⁠, the free-swimming Strouhal number decreases to $StA=0.28$⁠, because the normalized viscous shear drag (i.e., viscous drag coefficient) on the body is reduced, and the free-swimming velocity is increased for the same tail-beat amplitude. This also reduces the magnitude of the lateral momentum imparted in the wake as compared to the axial momentum and resulting in the wake at $ReL=50,000$ that exhibits a narrower spreading in the lateral direction than for the case with $ReL=5000$⁠. These simulations therefore highlight the role that size (or scale) has on the hydrodynamics of carangiform propulsion. The spreading angle of the wake has implications for the hydrodynamic interactions in fish schools, which are the topic of Sec. 3.1.3. We note here that the high Reynolds number simulations in this section as well as the later Sec. 3.3.2 on external aerodynamics indicate that IB methods, especially of the sharp-interface variety, are not limited to low-Reynolds number flows. Indeed, in a previous review article, Mittal and Seo [9] addressed the application of IB methods to high-Reynolds number flows as well as turbulent and high-speed flows. The rapid advancements in graphical-processing unit based acceleration of IB solver [56,57] are allowing these methods to further extend their reach into the higher Reynolds number regime, which has typically been outside the range of these methods.

The collective behavior of fish is a well-documented behavior observed in various species, characterized by the coordinated movement of fish in groups [58–60]. This collective behavior provides several advantages [61–66], including enhanced foraging efficiency, improved predator protection, stealthy swimming, and increased hydrodynamic performance. To investigate the dynamics of fish schooling and quantify the associated hydrodynamic and acoustic benefits, 3D computational simulations offer a robust alternative to traditional experimental approaches. The limitations of live animal studies, particularly due to the uncontrollable nature of real fish schools, highlight the necessity for simulation-based investigations [67]. In this context, the ghost-cell based sharp-interface immersed boundary method (vicar3d) employed here, that can accurately replicate the complex movements of fish, provides a powerful tool for modeling the hydrodynamic characteristics of swimming fish.

To investigate the benefits of fish schooling, we construct a 3D model of a fish based on the common mackerel (Scomber scombrus), consisting of a solid body and a membrane tail (Fig. 10(a)). The tailbeat dynamics are governed by well-established kinematic modes [68], ensuring a realistic simulation of the swimming behavior (Fig. 10(b)). The configuration and parameters of the swimming fish are based on the previous work of Seo and Mittal [55]. The fish is tethered in place in an incoming flow, $U=0.48Lf$⁠, and the Reynolds number based on the fish length and tail-beat frequency (⁠$L2f/ν$⁠) was 5000. This choice of parameters ensures that the total drag and thrust on the fish are balanced out and the fish is swimming at its terminal velocity.

Starting with simulating the solitary swimming fish (Fig. 10(c)), we extend our simulations to model schools of fish with varying configurations, focusing on the effects of relative positions and tailbeat phases. Contrary to the intuitive expectation that a school of fish would generate more chaotic and less efficient flow fields, our simulations reveal that a school could achieve significant hydrodynamic advantages with an appropriate configuration. Specifically, fish configured in optimal relative phases can utilize the leading-edge vortices generated by upstream ones, resulting in higher swimming efficiency [55,69]. For instance, a nine-fish school, shown in Fig. 10(e), could result in about 10% higher thrust per fish compared to the solitary fish swimming. Furthermore, our simulations demonstrate that coordinated tailbeat phases can reduce acoustic signatures, effectively allowing the school to swim more stealthily [65]. Figure 10(d) shows the optimal schooling configuration for a four-fish school to reach an about twofold deduction in the per-capita swimming noise. These findings emphasize the importance of precise spatial and temporal coordination in maximizing the benefits of fish schooling and demonstrate the robustness of the present ghost-cell sharp-interface immersed boundary method vicar3d [21].

In this section, we highlight several applications of the sharp-interface IB solver vicar3d [21] to physiological flows.

A thromboembolic event associated with left ventricular thrombosis (LVT) is one of the most critical complications for patients with myocardial infarction [70], and left ventricular (LV) flow patterns have been shown to play an important role in the formation of LVT. After a myocardial infarction, abnormal LV motions, LV remodeling, and the changes in ejection fraction combine in a complex way to affect LV flow patterns and increase the risk of LVT [70]. The high degree of nonlinearity associated with the flow dynamics makes it difficult to find direct correlation between the above LV abnormalities and the resulting LVT risk. Therefore, a direct assessment of LV flow is necessary for the identification of abnormal flow patterns or metrics associated with the LV flow that could serve as reliable indicators of the risk of LVT [71].

Based on the sharp-interface immersed boundary method framework of vicar3d, we developed a coupled chemofluidic computational model to investigate the quantitative correlation between LV flow and thrombus formation in modeled LVs. LV blood flow is simulated by solving the incompressible Navier–Stokes equations using the sharp-interface IBM, and this is coupled with biochemical modeling of thrombus formation which includes the coagulation cascade, fibrin polymerization, and platelet activation [72]. The method has been applied to canonical as well as patient-specific cases to investigate the correlation between flow metrics and LVT risk. Medical imaging data from LVT patients were obtained from the database at Johns Hopkins Hospital. Time dependent, 3D computational models of patient specific LVs are generated by a multimodality image-based reconstruction process.

The results of coupled chemofluidic simulations for four patient-specific cases are presented in Fig. 11 where the time-averaged flow patterns and the spatial distributions of thrombin are shown. One can see that the size and shape of the LVs are very different from case to case. Each case shows distinct hemodynamic flow pattern, which is caused by the combination of LV shape and motion. The simulation results show that if the blood flow from the mitral inlet does not propagate down to the apex, the incoming blood flow mostly resides in the basal region, and thus the flow strength near the apical infarction region is extremely weak. This results in poor apical washout and higher thrombin concentration, which leads to the higher LVT risk. Results from these simulations are described in more detail in Refs. [71] and [72].

While oral administration is the most common method of drug delivery to humans, it is a very complex process for an active pharmaceutical ingredient (API) to enter the body. This complexity is associated with the fact that drug absorption in the gastrointestinal tract is influenced not only by the drug and its formulation but also by the stomach's contents, motility, and the related fluid dynamics. The pressure and shear forces generated by stomach contractions, along with buoyancy effects, can lead to intricate pill trajectories, variable dissolution rates, and uneven emptying of the drug into the duodenum. To address these complexities, we developed a computational model using a sharp-interface immersed boundary, fluid–structure–chemical interaction solver, simulating drug dissolution in the physiological human stomach [73]. The model accurately represents the flow of gastric content due to stomach motility, resolves the six-degree-of-freedom motion of the drug in the gastric environment, and incorporates both diffusive and convective transport of the API concentration, which are critical to the drug dissolution process. Further details of the modeling procedure and the key model parameters can be found in Ref. [73].

Figure 12 shows the dissolution of the pill (surface erosion) simulated by directly solving the convection and diffusion equation for the API concentration. The volumetric distributions of dissolved API concentration are plotted at 10 s interval. In the simulation, the pill density is set to 1.2 times of the gastric fluid. The antral contraction wave of the stomach affects not only the motion of the pill but also the characteristic flow pattern called the retropulsive jet, which directly affects the convective transport of the API. Due to the higher density, the pill settles near the pylorus. The dissolved API concentration is then transported by the retropulsive jet to a region superior to the antral contraction wave and toward the body. The API is then mixed in the antrum by the retropulsive jet and the recirculating flow. The simulation results demonstrate that stomach motility and flow pattern have to be considered for the analysis of drug dissolution under physiological conditions. The developed model has also been applied to study the effects of body posture and gastroparesis on drug dissolution [74].

Flutter in the transonic regime is characterized by several nonlinear flow features including shock–boundary layer interactions, shock–vortex dynamics, and flow separation. These phenomena are known to contribute to a reduction in the critical flutter speed, referred to as the transonic dip [75]. The increased susceptibility toward wing flutter in this regime is a significant concern for the design and operation of military aircraft [76]. This is especially true due to the prominence of limit cycle oscillations in transonic flutter (opposed to flutter divergence), which can lead to fatigue damage to the aircraft structure and targeting issues. Transonic flutter can manifest in various forms including classical transonic buzz [77], shock-stall flutter [78], and boundary layer transition [79].

Simulating flutter in the transonic regime is especially challenging due to the requirement to accurately predict flow separation and shock locations, thus highlighting the need for high-order methods. Recently, we have successfully utilized vicas3d, a high-order sharp-interface IB solver [35], to investigate the instability mechanisms associated with single degree-of-freedom pitching flutter of a NACA0012 airfoil [36]. It was shown that the main energy extraction mechanism was flow separation induced by a propagating lambda-shock. Figure 13 shows a recent investigation of a three-dimensional NACA0012 airfoil undergoing forced vibration with pitching amplitude $Aθ=10 deg$ at normalized frequency $f*=0.1$⁠. We note that the use of forced vibrations as a means to study free vibration is well established in the arena of fluid–structure interaction [80,81] and combined with the notion of “energy-maps” [82], has been highly insightful for examining aeroelastic flutter of wings. The freestream Reynolds number is $Re∞=10,000$ and Mach number $M∞=0.7$⁠. Periodic boundary conditions are implemented in the spanwise direction. A total of 50.7 × 10⁶ mesh cells are used for the simulation, which is run on 960 processor cores. Figure 13(a) shows the pitch angle and corresponding power (normalized) extracted from the flow by the forced vibration. Instances with positive power indicated energy transfer from the flow to the foil which promote flutter. The flow field is visualized by isosurfaces of Q-criterion in Fig. 13(b) at $t*=33.5$⁠.

Wake flows due to bluff bodies such as cylinders and spheres are one of the most fundamental and well-studied problems in the field of fluid mechanics with relevance in a wide variety of engineering applications. Nevertheless, there remains limited knowledge regarding the impact of compressibility on these flows, mainly due to the lack of application-driven incentives until recent times. Prominent examples include re-entry vehicles, aircraft operating in extraterrestrial environments, and multiphase particle flows, e.g., rocket exhaust plume–surface interactions and ice crystal impact on hypersonic vehicles. Although a significant effort has been made previously to quantify drag characteristics (e.g., Ref. [83]), key features such as transition, instability modes, and wake dynamics require further investigation for such flows.

Using vicas3d, we have studied the complex wake dynamics of both cylinders and spheres in the supersonic regime [35]. An example, demonstrating the capability of the method is shown in Fig. 14. Contours of normalized density gradient magnitude (⁠$|∇ρ*|$⁠) are shown for an array of five 2D cylinders at a freestream Reynolds number of $Re∞=200,000$ and Mach number $M∞=1.7$⁠. A total of $5000×3600=18$M cells are used for the simulation. A weighted least squares minimization approach is used to implement the immersed boundary condition with a third-order accurate interpolation polynomial [33]. A hyperviscosity approach is also implemented for shock capturing [84].

Drones have become increasingly popular and are now integral to various fields, including entertainment, environmental monitoring, package delivery, military operations, and even the transport of critical medical supplies to remote areas [85]. They offer a more environmentally friendly alternative for material transport compared to traditional methods. However, the noise they produce remains a significant barrier to their widespread adoption in many applications [86].

The noise generated from helicopter and drone rotors has been a persistent problem for several decades, leading to several rotor design modifications to mitigate this noise [87]. An interesting attempt at rotor noise reduction was by Vogeley [88] almost 75 years ago, which involved replacing the conventional propeller of an airplane with a propeller with a much larger wetted area. This propeller was able to generate a similar thrust to conventional propeller while rotating at a lower speed, and the reduced rotational velocity resulted in a lower aeroacoustic noise. We have used a similar propeller design to model our quadcopter rotors, and the simulation using vicar3d is shown in Fig. 15 with Reynolds number based on tip velocity and radius of the rotor set to 4050 and each adjacent propeller rotating in the opposite direction to cancel out the net torque.

We have used the Ffowcs Williams and Hawkings equation [89] to calculate the loading noise from this quadcopter. The quadruple source term is ignored because of low surface Mach number, and the thickness noise is also ignored since we have used a zero thickness control surface for the rotors. The noise is dimensionalized for a rotor radius of 6 in and a rotation frequency of 40 Hz and is calculated at a distance of 5 m from the center of the quadcopter. The directivity plot is shown in Fig. 16, and it demonstrates the ability of the sharp-interface IB method to enable the exploration of the aeroacoustics of a full multirotor drone configuration.

The paper explores the developments and key features of sharp-interface immersed boundary methods, which have seen significant advancements over the past 25 years. These methods achieve a level of accuracy and precision in applying boundary conditions comparable to that of equivalent body-conformable grid methods. The paper outlines five criteria that must be met for an IBM to be considered a sharp-interface method. It also delves into two sharp-interface IB methods that have been the focus of our research group: the cut-cell method and ghost-cell based methods. The versatility of these methods is demonstrated through a wide array of flow simulations, ranging from biolocomotion and physiology to compressible aerodynamics and aeroacoustics. This highlights the broad applicability and effectiveness of these advanced immersed boundary methods.

Computational resources for this work were provided by the Advanced Research Computing at Hopkins (ARCH) core facility and by the National Science Foundation. We also acknowledge Dr. Dan Riskin and Dr. Sharon Swartz from Brown University for providing experimental data on bat wing kinematics.

• National Science Foundation (Grant Nos. CBET-2034983, CBET-1738918, PHY-1806689, DUE-1734744, CBET-2019405, CBET-2011619, CBET-1357819, and OAC1920103; Funder ID: 10.13039/100000001).

• Air Force Office of Scientific Research (Grant Nos. FA9550-21-1-0286 and FA9550-23-1-0010; Funder ID: 10.13039/100000181).

• Office of Naval Research (Grant Nos. N00014-22-1-2655 and N00014-22-1-2770; Funder ID: 10.13039/100000006).

• Army Research Office (Cooperative Agreement No. W911NF2120087; Funder ID: 10.13039/100000183).

• National Institutes of Health (Grant No. 5R21GM139073-02; Funder ID: 10.13039/100000002).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.