The flow field is modeled using unsteady incompressible Navier-Stokes equations: ∂ u i ∂ x i = 0 ∂ u i ∂ t + ∂ u i u j ∂ x j = − ∂ p ∂ x i + 1 R e ∂ 2 u i ∂ x j ∂ x j where u i represent the velocity components in the three directions and p is pressure. R e is the Reynolds number, defined as U ∞ ν , where U ∞ is the incoming flow speed, D is the characteristic length, and ν is the kinetic viscosity of water at 20°C. The presence of the whiskers in the flow field is ignored in this study.

The current study employs an in-house immersed boundary method based incompressible flow solver, as described by Mittal et al. (2008). The solver discretizes the spatial terms using a cell-centered collocated arrangement of the primitive variables u i and p . To integrate the equations in time, a three sub-steps fractional step method based on Van Kan (1986) is employed. In the first sub-step, an intermediate velocity field, u * , is determined by solving a modified momentum equation without the pressure term. The convective terms in this equation are discretized using a second-order Adam-Bashforth scheme, while the diffusion terms are discretized using an implicit Crank-Nicholson scheme to remove the viscous stability constraints. In the second sub-step, the pressure corrections are computed by solving a pressure Poisson equation, and in the final step, the pressure field and velocity field are updated using the obtained pressure corrections. To treat the immersed boundaries, a multidirectional ghost-cell methodology is employed, resulting in both local and global second-order accuracy. The solver used in this study has been successfully applied to various biological applications, including human phonation, insect flight, fish swimming, and seal whiskers inspired flow sensing (Geng et al., 2017; Liu et al., 2018; Liu et al., 2019b; Bodaghi et al., 2021a; Bodaghi et al., 2021b; Liu et al., 2021). Details of the solver can be found in Mittal et al. (2008).

Each whisker is modeled as elastic structure and dynamics of whisker is governed by the Navier equation: ∂ 2 d i ∂ t 2 = 1 ρ w ∂ σ i j ∂ x j + f i where d i , ρ w , σ i j and f i are displacement vector, whisker density, stress tensor and body force vector, respectively. An in-house finite element method code is utilized to solve the equation for d i based on the hydrodynamic loads obtained from the flow solver. This solver was fully validated by Liu et al. (2019a) and was successfully employed for different biological simulations, including fish swimming (Liu et al., 2018) and seal whiskers (Liu et al., 2021).

Upon solving the dynamic equation at each time step, the hydrodynamic loadings along each whisker are estimated based on the local flow speeds using a simplified drag-lift model. To achieve this, each individual whisker is discretized into a series of whisker segments, each one spanning one segment. The hydrodynamic force acting on each segment is computed using F D , F L = C D , C L ∙ 0.5 ρ U 2 A r e f , where F D and F L are the drag and lateral forces along the primary and minor axis of the cross section, respectively. ρ is flow density, U is local flow speed and A r e f is referenced cross section area. C D and C L are the drag and lateral force coefficients of the whisker segment and depend on local Reynolds number and angle of attack of the segment. We obtained these coefficients by using a tabulated function from our previous numerical parametric study. The resulting hydrodynamic loadings are then applied the centroid points of each segment, together driving the motion of the entire whisker. The details of hydrodynamic load calculation can be found in Liu et al. (2021).

The seal whisker array model, depicted in Figure 1A, is adopted from our previous study (Liu et al., 2021), which closely replicates a real harbor seal’s whisker array. The model comprises of two arrays, each containing 44 individual whiskers, placed on the muzzle of a realistic harbor seal head model. The length, curvature, tapering, location, and orientation of each whisker are taken from the reported data of realistic harbor seal whisker arrays (Murphy, 2013; Graff, 2016) and are illustrated in Figures 1B, C. The undulated geometry of each whisker is constructed using a simplified elliptical cylinder model, superimposed with two sinusoidal undulations with a 180° phase shift of the axes of the elliptical cross-sections (Figure 1D) (Rinehart et al., 2017). The material properties of the whiskers are also determined using experimental data (Hans et al., 2014). For more detailed information, please refer to Liu et al. (2021).

Figure 2 illustrates the simulation setup. A circular plate with a diameter of 40 mm is placed 30 cm in front of the nasal tip of a seal head. This diameter is determined based on previous studies (Wieskotten et al., 2011; Liu et al., 2021) for generating typical-size vortices of prey. The seal head and the plate are immersed into a 900   m m × 1000   m m × 950   m m flow domain. A uniform flow speed of 20 cm/s from left to right is applied as far field conditions on left, top, bottom, front and back boundaries of the flow domain and the outflow condition is implemented on the right boundary of flow domain with zero velocity gradient and zero pressure gradient. To reduce the computational time, we increased the kinematic viscosity of fluid to 4 × 10 − 5 m 2 / s , leading to a reduced Reynolds number of R e D = 200 , defined using the diameter of the circular disk. To ensure grid independence, two different grids were tested, with resolutions of 128 × 128 × 256 and 256 × 256 × 256 grid cells in the x, y, and z directions, respectively. The highest grid density is concentrated near the circular disk, its wake region as well as the whisker arrays, with minimum grid cell sizes of 3.2   m m × 4.4   m m × 2.6   m m and 1.6   m m × 2.2   m m × 1.3   m m for the coarse and fine grids, respectively. The grid independence study shows that using the coarse grid results in less than 0.5% difference in the drag force experienced by the circular disk and 1% difference in the vortex shedding period compared to the fine grid. Considering the significant reduction of computational cost, the 128 × 128 × 256 coarse grid is chosen for the study, as depicted in Figure 2B. The flow simulations are executed on the XSEDE Expanse cluster (AMD EPYC 7742 type CPU, clock speed: 2.25 GHz, and flop speed: 4608 GFlop/s) using 16 processors. The simulations are run for 5,000 time steps to capture three steady vortex shedding cycles, resulting in an average computational cost of 224 CPU hours per case for the flow simulation.

The whiskers are discretized using ten-node tetrahedral cells, resulting in a range of 9,700–34,100 cells depending on the length of the whiskers. The whiskers are fixed at the base, preventing any movement or rotation. The solid simulations are performed on the same XSEDE Expanse cluster, with one processor allocated for each whisker. Thus, for each parametric case that includes 88 whiskers, the average computational cost for the solid simulation is 184 CPU hours.

A group of parametric simulations were conducted by systematically varying the location (X, Y) and orientation (α, β) of the circular plate, as defined in Figure 2A. The ranges of these parameters were determined to ensure the wakes would intersect at least one of the whiskers. The selected parameters are: X from −120 mm to 120 mm with an interval of 60 mm, Y from 0 to 200 mm with an interval of 50 mm, with an interval of 50 mm, α from −45° to 45° with an interval of 22.5°, and β from −45° to 45° with an interval of 22.5°. This results in a total of 625 simulation cases.

Vision Transformer (ViT) (Vaswani et al., 2017; Dosovitskiy et al., 2020) has demonstrated effective structural interpretation of spatial correlations between pixels in images by dividing raw images into multiple image patches and leveraging comprehensive imaging information. Flexibility of ViT to handle various input sequences allows for the utilization of video information in addition to image information. In this study, we propose a network architecture called Video Transformer (VT) as a straightforward extension of ViT, designed specifically for video processing. VT initially splits the video into spatio-temporal tokens and applies a self-attention model using spatial and temporal blocks to extract features based on spatial and temporal changes of moments. These features are then optimized by passing through transformer blocks to determine the target position and orientation. In this section, we first introduce ViT, followed by a discussion on tubelet embedding for video processing (Arnab et al., 2021), and finally, we present the self-attention model for video processing (Arnab et al., 2021).

The primary input signal for the network is derived from the time history of the bending moment at the base of each whisker. The total bending moment is decomposed in two components: M x ′ and M y ′ , where x ′ and y ′ denotes body-fixed coordinates of each whisker, as depicted in Figure 1C. Based on the direction of x ′ and y ′ , M x ′ represents the bending moments resulting from the lift force while M y ′ represents the bending moments resulting from the drag force. Please note that an additional simulation of the freestream passing the seal head without an upstream plate was conducted and the bending moments obtained from this simulation were then subtracted from the bending moments observed in the other cases where wakes were present. The process is to eliminate the influence of the freestream and solely examine the effect of the wakes.

Tactile vibrissal systems in general possess two distinct types of mechanoreceptors in their whiskers, namely Merkel cells (SAI mechanoreceptors) and Pacinian corpuscles (PC mechanoreceptors), specialized in detecting static and dynamic changes (Zimmerman et al., 2014). Merkel cells are adept at sensing static stimuli, such as sustained pressure, while Pacinian corpuscles excel at perceiving dynamic changes, including vibrations and rapid movements. To emulate the underlying nerve system, we further decompose the bending moments into two components: time-averaged (DC) components, denoted as M D C , and oscillatory components (AC), denotes as M A C . This yields a four-component vector: ( M D C , x ′ ; M A C , x ′ ; M D C , y ′ ; M A C , y ′ ), obtained from each whisker for the network input.

The subsequent step involves a data transformation process aimed at consolidating the inputs from individual whiskers to generate a signal map, as depicted in Figure 4. Initially, the base locations of each whisker were extracted from the model and then topologically mapped onto a 9 × 22 two-dimensional rectangular grid based on their relative positions. Subsequently, each grid point was treated as a pixel and associated with the bending moment vector of the corresponding whisker. Each moment component was considered a distinct channel, with the intensity of each channel determined by the magnitude of its respective moment component. As a result, a “representative sensing image” comprising four channels was created. In cases where a pixel lacked correspondence to any whisker, a NULL (0) value was assigned. Given that bending moments exhibit temporal variations, the resulting data ultimately takes the form of a video file, reflecting the dynamic nature of the sensory input.

The network training process involves regression training for 1,000 epochs using the Adam optimizer. A batch size of 25 and an initial learning rate of 0.001 are used, along with a multi-step learning rate scheduler. The input samples are divided into train and validation groups using the 5-fold method, with 500 samples for training and 125 samples for validation. Each video sample has a size of 80 × 9 × 22 × 4 , representing the temporal dimension (80, corresponding to the last 1,500 flow time steps to only consider the steady vortices), the whisker array dimensions ( 9 × 22 ), and the moment channels (4). To capture the relationships between each whisker, the model parameters are initialized as follows: t ′ = 8 , h = 1 and w = 1, with L = 4 and G = 96.

The implementation of the method utilizes PyTorch 1.12.1 and Skorch 0.11.0 (Tietz et al., 2017) frameworks. The training process is performed using an NVIDIA GeForce RTX 2080 Ti GPU with 11GB memory.

Figure 5 shows the wake structure induced by the plate in the exemplary cases, determined by the λ 2 criterion, viewed from different viewpoints: 3D, top view, and side view. To aid the comparison, the baseline case is added in the figure. Except for α 45 β 22.5 case, the wake displays two distinct vortex rings that alternate between the upper and lower sides ( X 60 , Y 50 , and β 22.5 cases) or the left and right sides ( α 22.5 case), with a phase difference of 180°. In the baseline case, where the plate is directly in front of the seal’s nose, the wake interacts with all whiskers. As the plate moves upward in the X 60 case and rightward in the Y 50 case, the wake translates accordingly without changing the pattern. In the α 22.5 case, where the plate tilts to the left, the wake inclines in the same direction. Additionally, due to plate rotation, the vortex ring forming on the left side is stronger than the one on the right side. In the β 22.5 case, where the plate tilts upward, the wake also inclines upward accordingly, and the rotation generates a stronger vortex ring on the upper side than the lower side. In contrast, the α 45 β 22.5 case does not exhibit vortex shedding but instead displays two streamwise vortex legs attached to the plate’s sides.

The vortex shedding frequency also changes with the orientation of the plate. To illustrate the relationship, Figure 6A presents the vortex shedding frequency in the parametric space of α and β. When vortex shedding is absent, a zero value is assigned. The plot reveals that the vortex shedding frequency increases symmetrically along both the α and β axes. When one orientation angle is at its maximum while the other is non-zero, no vortex shedding occurs. To account for the changes in both the vortex shedding frequency and the frontal area of the plate, the Strouhal number was calculated for each case using the formula S t L ′ = f L ′ / U ∞ , where f represents the vortex shedding frequency, U ∞ denotes the inflow velocity, and L ′ is the hydraulic diameter of the frontal area of the plate. Figure 6B demonstrates the change in the Strouhal number with orientation, mirroring the behavior of the vortex shedding frequency. The maximum increase in the Strouhal number is 80%. A similar observation was reported by Chen and Fang (1996), where a decrease in the angle of attack of an inclined plate (comparable to increasing α and/or β in our study) resulted in an increase in the Strouhal number of its vortex shedding.

We partition the whisker arrays into four regions: LT (left top), LB (left bottom), RT (right top), and RB (right bottom) whiskers, and calculated the average signals, M x ′ and M y ′ , across the whiskers within each region. Figure 7 presents the phase-averaged time history of M x ′ and M y ′ in each region in each exemplary case. Note that M x ′ is mostly negative values in all the cases, while M y ′ exhibits both positive and negative values. The behaviour of the signals aligns with the temporal and spatial dynamics of the vortices in each case, (e.g. Frequency). In the X 60 case, the signals are evenly distributed across the regions since the vortices impact the seal head at the center. In the Y 50 case, where the wake translates rightward, both M x ′ and M y ′ are stronger on the right whiskers. The right bottom whiskers experience stronger signals compared to the right top whiskers because of the impact location of the vortices. In the α 22.5 case, the left whiskers experience stronger signals than the right whiskers, consistent with the plate rotation that generates stronger vortices on the left side. Additionally, the left bottom whiskers experience stronger moments than the left top whiskers. In the β 22.5 case, interestingly, despite the upper vortex ring being stronger than the bottom one, the bottom whiskers experience stronger M x ′ than the top whiskers. This is attributed to the upward inclination of the wake, causing the upper vortex ring moving away from the whiskers. Also interestingly, different from M x ′ , M y ′ shows no variation across the regions in this case. In the α 45 β 22.5 cases, both M x ′ and M y ′ are almost constant with time while the left bottom whiskers have stronger moments compared to the other whiskers due to the stronger wake in this region.

Overall, the signals at the root of the whisker arrays exhibit temporal and spatial patterns that align with the dynamics of the vortices. There is a clear correlation between the arrival of the vortices and the temporal variation in the signals. The regions where the vortices directly impact the whiskers tend to exhibit stronger signal amplitudes. The distinct vortex dynamics in each case result in different temporal and spatial patterns of the signals, providing an explanation for the discernible differences among the cases.

The neural network was trained successfully to predict the location and orientation of the upstream plate. Figure 8A illustrates the best mean squared error (MSE) losses of the training and validation groups over 1,000 epochs, indicating the convergence of the network. The final MSE loss for the training and validation corresponding to the best fold is 0.0255 and 0.1064, respectively. The decision to limit the epochs to 1,000 was based on the observation that the drop in MSE loss for both groups during the last 100 epochs is less than 0.1%. The accuracy of the trained network corresponding to the best fold is depicted in Figure 8B through a comparison of the ground truth and predicted locations and orientations of the upstream plate in the validation cases. The accuracy metric is quantified using the following formula: A c c u a r c y = 1 − ∑ y G T − y Pr R a n g e N × 100 where y G T represents the ground truth value, y Pr denotes the predicted value, R a n g e represents the range of changes of y G T , and N is the number of validation cases. The obtained 5-fold-average accuracy of the parameters X , Y , α , and β is 91.4 ± 1.1 % , 91.0 ± 1.1 % , 84.8 ± 2.6 % , and 83.3 ± 0.6 % , respectively, resulting in an average accuracy of 87.6 ± 0.7 % for the network.

An essential capability of the interpretable network is its ability to generate temporal and spatial importance of the signals within the trained network, allowing for the identification of significant temporal-spatial signal patterns that contribute to the network’s predictions. In our study, as outlined in Section 3.3, we utilized a temporal dimension of 80 and a time interval of t = 8, resulting in 10 distinct time patches for the prediction. In the X 60 and Y 50 cases, these 10 time patches correspond to two consecutive cycles, while in the α 22.5 and β 22.5 cases, they correspond to three consecutive cycles.

In Figure 9, we present the temporal importance (TI) analysis of the signals in each representative case. Except for the α 45 β 22.5 case, the results highlight that the temporal importance is not uniformly distributed across all time instances. Instead, certain time points exhibit significantly higher importance values, while others demonstrate relatively lower importance values. This disparity suggests that specific time instances play a crucial role in the prediction process, while others have a lesser impact. Furthermore, it is noteworthy that the temporal importance patterns differ across all the cases. Each case exhibits a unique distribution of importance values across the time patches. In the α 45 β 22.5 case where vortex shedding is absent, the temporal importance is nearly uniform acorss the time patches within a cycle. It is because the flow field remains nearly unchanged with time.

The variation in temporal importance patterns highlights the importance of considering the temporal dynamics and their specific correlations with the underlying flow structures when interpreting the prediction process. To illustrate this relationship, we present the vortex structures determined by the λ 2 criterion and contours of flow velocity magnitude at two planes intersecting the whisker arrays at the first five individual time patches in each representative case in Figure 10. Note that the flow velocity in the case of freestream passing the seal head without an upstream plate was subtracted from the flow velocity in other cases where wakes were present to focus on the disturbance induced by the wake. To facilitate analysis, the temporal importance map is also included in each case. These five patches cover either one or more than one cycle of the wake. The α 45 β 22.5 case is not included as the temporal importance and flow structure remains unchanged in this case.

In the X 60 case, time patch 1 exhibits significantly higher importance compared to the other time patches. This observation corresponds to the vortex dynamics, where at time patch 1, the lower vortex ring actively interacts with the whiskers. In contrast, at the other time patches, either the top vortex ring arrives in the whisker area but lacks strong interaction due to its distance from the whiskers, or the bottom vortex ring has already moved past the whiskers. The flow velocity contour reveals that the flow disturbance is strongest at time patch 1 around the whisker compared to other time patches.

In the Y 50 case, time patch 4 demonstrates significantly higher importance compared to other time patches. This time patch corresponds to the instances when the top vortex ring interacts with the right whiskers. At the other time patches, either the top vortex ring has not yet reached the whiskers or has already moved past them. In this case, the bottom vortex ring has minimal effect due to its distance from the whiskers resulting from the translation of the upstream plate. The flow velocity contour also shows that the vortex ring interacts with the right whisker array and the disturbance is the strongest on most whiskers at time patch 4.

In the case of α 22.5 time patch 4 exhibits the highest importance, while time patches 1 and 5 are slightly less important, and time patches 2 and 3 have minimal importance. Interestingly, upon analyzing the vortex dynamics, it is observed that time patches 2 and 3 capture the interactions between the left whiskers and the stronger vortex ring from the left side. Conversely, time patch 4 captures the interactions between the right whiskers and the weaker vortex ring shed from the right edge. Time patches 1 and 5 represent the transitional stages, where the weaker vortex ring is about to leave and the stronger vortex ring is about to arrive. The flow velocity contour further confirms that the flow disturbance is strongest at time patches 2 and 3 due to the presence of the stronger vortex, whereas it is weaker at time patch 4 attributed to the weaker vortex ring. One possible explanation for the higher importance of the weaker vortex in this case is that the stronger vortex ring shed from the left edge aligns towards the center of the nose, similar to the vortices observed in other cases. In contrast, the weaker vortex ring is oriented perpendicular to it and only interacts with the right whisker arrays. This unique interaction provides additional crucial orientation information, making the weaker vortex more influential in the prediction process.

In the case of β 22.5 , time patches 1, 4, and 5 exhibit greater significance compared to time patches 2 and 3. Through vortex dynamics analysis, it is observed that time patches 1 and 5 capture the arrival of the top vortex at the whiskers, while time patch 4 captures the interactions between the bottom vortex and the whiskers. Interestingly, time patch 2 captures the strong interactions between the top vortex ring and the whiskers, but its contribution is nearly negligible. This observation is further supported by the flow velocity contour. The flow disturbance during time patches 1, 2, and 5 is primarily associated with the top vortex, whereas during time patch 4, it is linked to the bottom vortex. Time patch 3 presents a scenario where disturbances from both vortex rings are present. Notably, similar to the α 22.5 case, the two vortex rings differ in strength, with the top ring being stronger than the bottom one, but the weaker ring plays the most important role in prediction. It is possible for the same reason that the weaker ring possesses a unique orientation and only interacts with the bottom whiskers. However, unlike the α 22.5 case, the significance of the other ring is increased in this case. It may be because the orientation information provided by the bottom ring is insufficient for accurate prediction, thereby necessitating the incorporation of additional information.

The current interpretable deep learning network also provides the importance of each whisker through spatial importance (SI) maps at individual time patch. Figure 11 illustrates the locations of the whiskers with the highest importance values during each of the first five time patches in each case. These whisker locations generally exhibit a strong correlation with areas of strong flow disturbance, as depicted in Figure 10. Furthermore, the spatial distribution of these important whiskers provides insights into the location and size of the vortex rings in the wake. For instance, in the X 60 case, the important whiskers reflect the vortex ring impacting the center of the seal head. In the Y 50 case, the important whiskers correlate with strong flow disturbance on the right side, and their extension reflects the extension of the vortex ring. In the α 22.5 case, the important whiskers capture both vortex rings, with one located on the left side and the other on the right side. Also, they reveal that the extension of the left vortex ring is larger than the right one. Similarly, in the β 22.5 case, the important whiskers capture two vortex rings. In this scenario, they capture that the higher vortex ring is larger and situated near the center, while the lower vortex ring is smaller and positioned closer to the bottom.

To gain a deeper understanding of the signals used for prediction, Figures 12, 13 present the maps of M x ′ and M y ′ on the whisker arrays during each time patch for each case. Note that M x ′ represents the moment generated by the lift force, while M y ′ represents the moment generated by the drag force. A strong correlation between the important whiskers and M x ′ maps is observed across all cases, indicating the significant role of M x ′ in the prediction process. Conversely, we did not observe any notable correlation with the M y ′ maps, suggesting that M y ′ has minimal impact on the prediction. Interestingly, the direction of M x ′ plays a crucial role in the prediction, as the important whiskers consistently correlate with positive M x ′ values across all cases or, in cases where positive M x ′ values are absent, with less negative M x ′ values. It is important to note that in our algorithm, the calculation of M x ′ involves subtracting the moment induced by the freestream flow. This subtraction is intended to focus on the moment induced specifically by vortices. Further studies are needed to clarify whether the direction information remains important if the original moment, without subtracting the freestream moment, is used. Nevertheless, the results suggest that the bending moment generated by the lift force plays a crucial role in the prediction process, and the direction of this moment may also be important.