JEB desktop wallpaper calendar 2016

Journal of Experimental Biology partnership with Dryad

Using computational fluid dynamics to calculate the stimulus to the lateral line of a fish in still water
Mark A. Rapo, Houshuo Jiang, Mark A. Grosenbaugh, Sheryl Coombs


This paper presents the first computational fluid dynamics (CFD) simulations of viscous flow due to a small sphere vibrating near a fish, a configuration that is frequently used for experiments on dipole source localization by the lateral line. Both two-dimensional (2-D) and three-dimensional (3-D) meshes were constructed, reproducing a previously published account of a mottled sculpin approaching an artificial prey. Both the fish-body geometry and the sphere vibration were explicitly included in the simulations. For comparison purposes, calculations using potential flow theory (PFT) of a 3-D dipole without a fish body being present were also performed. Comparisons between the 2-D and 3-D CFD simulations showed that the 2-D calculations did not accurately represent the 3-D flow and therefore did not produce realistic results. The 3-D CFD simulations showed that the presence of the fish body perturbed the dipole source pressure field near the fish body, an effect that was obviously absent in the PFT calculations of the dipole alone. In spite of this discrepancy, the pressure-gradient patterns to the lateral line system calculated from 3-D CFD simulations and PFT were similar. Conversely, the velocity field, which acted on the superficial neuromasts (SNs), was altered by the oscillatory boundary layer that formed at the fish's skin due to the flow produced by the vibrating sphere (accounted for in CFD but not PFT). An analytical solution of an oscillatory boundary layer above a flat plate, which was validated with CFD, was used to represent the flow near the fish's skin and to calculate the detection thresholds of the SNs in terms of flow velocity and strain rate. These calculations show that the boundary layer effects can be important, especially when the height of the cupula is less than the oscillatory boundary layer's Stokes viscous length scale.


All fishes possess a mechanosensory lateral line system, which responds to the surrounding water motion relative to the fish's skin. Such water motion can be generated due to a variety of biotic and abiotic events, including encounters with prey, predators, conspecifics and inanimate obstacles. Consequently, the lateral line system plays an important role in mediating fish behavior in different ecological contexts. Over the years, there is a rich body of literature documenting various aspects of lateral line function, ranging from biomechanical and neural principles of operation to behavioral significance (for reviews, see Dijkgraaf 1963; Kalmijn, 1988; Denton and Gray, 1988; Coombs et al., 1988; Coombs et al., 1992; Bleckmann, 1993; Montgomery et al., 1995; Schellart and Wubbels, 1997; Coombs and Montgomery, 1999; Bleckmann et al., 2001; Janssen, 2004; Mogdans, 2005; van Netten, 2006; Bleckmann, 2008).

The lateral line system consists of superficial (SN) and canal (CN) neuromast subsystems, which are morphologically, physiologically and functionally different (e.g. Coombs et al., 1988; Münz, 1989). In terms of the basic structure, both types of neuromasts consist of mechanosensory hair cells covered by a gelatinous cupula. However, SNs are generally smaller in diameter (≤100 μm) than CNs (up to several thousands of microns) and contain fewer hair cells (<∼100) compared with CNs (∼200–10,000) (Münz, 1989). Although the relative numbers and spatial distribution of SNs versus CNs over the head and body are quite species-specific, most if not all fish species (e.g. the mottled sculpin Cottus bairdi) have both types of neuromasts (Coombs et al., 1988). A reported exception is the plainfin midshipman fish, Porichthys notatus, whose trunk lateral line only has SNs (Weeg and Bass, 2002). The clear distinction between SNs and CNs is that SNs are located superficially on the skin surface whereas CNs reside just below the skin surface in a fluid-filled canal that opens up to the surrounding water via a series of pores – one pore between each pair of neuromasts. Thus, whereas SN cupulae protrude directly into the water surrounding the fish, CN cupulae protrude into the canal fluids. Hair cells of both types are mechanically excited by viscous coupling with the surrounding fluid motion via the overlying cupula. An SN cupula, protruding into the surrounding water, responds to the local fluid velocity field. A CN cupula, being inside the canal, responds to the fluid motion induced by the local net acceleration of the water against the fish skin, which is also proportional to the pressure gradient across the surrounding two canal pores (Denton and Gray, 1983; Kalmijn, 1989). Physiologically, SNs and CNs are likely to be innervated by different afferent fibers (Münz, 1985). Thus, fish have separate channels for measuring the velocity and acceleration fields produced by not only periodic (Denton and Gray, 1983; Münz, 1985; Coombs and Janssen, 1989; Coombs and Janssen, 1990; Kroese and Schellart, 1992) but any stimuli (Bleckmann, 2008).

Attention has been paid to the relative contributions of SNs and CNs to the orienting behavior of fish in response to both live and artificial prey (e.g. a chemically inert vibrating sphere). For the mottled sculpin, both neurophysiological and controlled behavioral studies suggest that the acceleration-responsive CNs, rather than the velocity-responsive SNs, mediate the orienting behavior (e.g. Coombs and Janssen, 1990; Coombs et al., 2001). However, sculpin larvae can feed in the dark on free-swimming Artemia at distance with the aid of SNs (Jones and Janssen, 1992). Despite the fact that the SNs of the sculpin larvae are eventually embedded into the canal to become CNs, the prey detection at the larval stage is due to movement of the SN cupulae in response to local flow velocities. Similar results were also found for the larvae of the willow shiner Gnathopogon elongatus caerulescens for whom it was shown that the number of Artemia consumed per larva was proportional to the cupular length and saturated at lengths above ∼100 μm (Mukai et al., 1994). Later, Mukai performed a controlled experiment using willow shiner larvae that demonstrated the role of mechanoreception by the SNs on prey detection (Mukai, 2006). Abdel-Latif et al. showed that, under still water conditions, the blind cave fish Astyanax mexicanus could still detect and approach a small vibrating sphere, presumably by means of its SN subsystem after the destruction of its CN subsystem (Abdel-Latif et al., 1990). However, there is a debate (e.g. Coombs et al., 2001; Mogdans, 2005) concerning the exact detection mechanism, as the pressure field around a dipole source could theoretically be detected by the pressure-sensitive ear of this species.

The stimulus to the lateral line system is the relative motion between fish skin and adjacent water, which is described by the spatio–temporal varying flow velocity and pressure fields surrounding the lateral line system. In general, such flow velocity and pressure fields are not (fully) measured in neurophysiological and behavioral studies – particularly in regions close to the skin where boundary layers are important. The inability to measure and specify the stimulus hinders our understanding of the lateral line system (Coombs and Montgomery, 1999). Potential flow theory (PFT) has been used in the past to address this issue. For example, the potential dipole source flow equations were used to model the pressure field due to a vibrating sphere near a fish body (e.g. Coombs et al., 1996; Coombs and Conley, 1997a; Coombs and Conley, 1997b; Conley and Coombs, 1998; Coombs et al., 2000; Curcic-Blake and van Netten, 2006). The same equations were used to calculate the slip flow velocity along the fish skin, caused by a nearby vibrating sphere (e.g. Kroese et al., 1978). (In PFT, the boundary condition at a solid wall requires only that the fluid velocity normal to the wall be equal to the wall velocity. Thus, the fluid is allowed to `slip' past the wall surface in the wall tangential direction.) This slip flow velocity was taken as the stimulus to the SNs, with the assumption that the cupular lengths were long enough to penetrate the viscous boundary layer along the fish's skin. More sophisticated potential flow solutions have been developed to calculate the slip flow velocity distribution (and pressure distribution) over idealized fish body geometry (Hassan, 1985; Hassan, 1992a; Hassan, 1992b; Hassan, 1993). All of these studies ignore the fact that the real flow has to satisfy the no-slip boundary condition at the fish skin due to viscosity, i.e. zero relative velocity at the skin in both the normal and tangential directions [see pp. 140–143 in Panton (Panton, 1996)].

Computational fluid dynamics (CFD) solves the Navier–Stokes equations, which include the viscous terms. CFD simulations have been previously employed to investigate tadpole swimming (Liu et al., 1996; Liu et al., 1997), fish undulatory swimming (Carling et al., 1998; Kern and Koumoutsakos, 2006; Borazjani and Sotiropoulos, 2008), dorsal–tail fin interaction in swimming fish (Akhtar et al., 2007), oral cavity flow in ram suspension-feeding fish (Cheer et al., 2001), jet flow behind a modeled swimming squid (Jiang and Grosenbaugh, 2006) and drag forces acting on a simulated neuromast inside a fish lateral line trunk canal with the canal flow driven by a two-dimensional (2-D) vortex street outside the canal (Barbier and Humphrey, 2008). To our knowledge, no previous CFD simulations have studied the flow due to a vibrating sphere (i.e. a dipole source) or a swimming prey-like object near a three-dimensional (3-D) fish body and calculated directly the stimuli to the lateral line system. CFD is useful because it has the ability to consider realistic fish body geometry and realistic spatial arrangements of the fish body and the signal source and because it can simultaneously output both the pressure and flow velocity fields at the lateral line locations. In the present study, we carried out a CFD numerical investigation of the flow due to a vibrating sphere near a mottled sculpin in still water. We performed a series of simulations of a prey-tracking sequence of a mottled sculpin as it responded to an artificial prey (i.e. vibrating sphere). The approximate shape and orientation of the fish body were taken directly from a video recording of a real tracking sequence (Coombs and Conley, 1997a). The series of simulations show explicitly how the presence of the fish body perturbs the contour lines of constant pressure of the dipole field.

The flow due to the vibrating sphere must satisfy the no-slip boundary condition at the fish's skin. This produces an oscillatory boundary layer of sheared flow that affects the magnitude of the hydrodynamic signals detectable by the SNs. As important as this effect is, we were not able to simulate it directly for a 3-D fish-shaped body because of limited computational resources. Instead, we used CFD to calculate the 3-D oscillatory boundary layer flow along a flat plate produced by a nearby vibrating sphere. We then used this result to validate an analytical model of an oscillatory boundary layer flow. The analytical model was then applied to a number of previous neurophysiological and behavioral experiments to calculate detection thresholds that take into account viscous effects.

This paper considers only flow due to the motion of the vibrating sphere in calm water. A future paper will include the added effects of an ambient, unidirectional current.


Simulations of flow around a sculpin due to a nearby vibrating sphere

Our computational domain and meshes for the fish simulations were constructed according to the experimental designs of Coombs and Conley (Coombs and Conley, 1997a; Coombs and Conley, 1997b) and by using the geometry and mesh generation software, GAMBIT (Lebanon, New Hampshire, USA). The geometry consisted of a mottled sculpin-like body, sitting on a solid flat bottom and a nearby 3 mm-radius solid sphere executing a vibration. The sphere performed sinusoidal motions along an axis for a series of bursts of 500 ms on and 500 ms off so as to dismiss the effects of acoustic streaming, which took more than 500 ms to develop. The sculpin body was constructed to approximately match the dimensions of the real fish. The length of the main body of the fish was 8 cm, with a portion of the tail extending an additional 1 cm. The body was 2 cm in width at its thickest point, which is just in front of the pectoral fins. The dorsal fin extended about 0.5 cm above the fish in its front and about 1 cm above in its rear. The pectoral fins and tail fin were 0.5 mm thick whereas the dorsal fin was 1 mm thick. Both a 2-D setup and a 3-D setup were considered. The purpose was to compare the 2-D results with the 3-D results and to find out if the 2-D setup was accurate enough to reproduce the involved fluid physics. (The 2-D assumption is often adopted in studies but it is necessary to check its validity before use.)

In the 2-D setup (Fig. 1), the cross-sections of the sculpin body and of the sphere formed the inner boundaries of a rectangular computational domain (275×100 cm). The sculpin body was held stationary. The cross-section of the sphere was given a prescribed axis of motion along with a solid-body vibration with the following time history: Math(1) where U0 is velocity amplitude of the sphere vibration, t is time and ω=2πf where f is the vibration frequency in Hz. A no-slip boundary condition was prescribed at the sculpin body and at the perimeter of the sphere cross-section. Also prescribed was a zero pressure inlet boundary condition at the outer boundaries of the domain (to approximate an infinite domain). The area between the inner and the outer boundaries of the domain was discretized into triangular meshes. The consecutive nodal points on the sculpin cross-section were 2 mm apart, close to the spacing between two consecutive lateral line CNs in the trunk canal of a real sculpin (Coombs et al., 1988). The perimeter of the sphere cross-section was equally divided by 20 nodal points. The size of the triangular meshes was increased gradually from the inner boundaries to the outer boundaries of the domain. A dynamic mesh model was used to represent explicitly the vibrating motion of the sphere cross-section, i.e. the meshes surrounding the sphere cross-section deform as the sphere cross-section moves (Fig. 1B).

In the 3-D setup (Fig. 2), a box-shaped computational domain (60×40×10 cm) was considered. The surface of the 3-D sculpin body was discretized into triangular meshes with 2 mm edge lengths. The sphere surface was divided into triangular meshes with 1 mm edge lengths. The volume between the sculpin body surface, the sphere surface and the surface of the box-shaped computational domain was discretized into tetrahedral control volumes. A deforming mesh zone was placed around the sphere surface to accommodate the sphere motion. Suitable boundary conditions were prescribed similar to those in the 2-D setup.

In the present still-water case, flow was generated exclusively from the vibration of the sphere. The flow was assumed laminar, incompressible and Newtonian and was governed by the unsteady incompressible Navier–Stokes equations together with the continuity equation. Throughout this study, the fluid density, ρ, was 1.0×103 kg m–3 and the fluid kinematic viscosity, ν, was 1.0×10–6 m2 s–1. The governing equations with the above-described computational domains and simulation setups were solved by a commercially available, finite-volume code, FLUENT™ (v. 6.2.16, Lebanon, New Hampshire, USA). The third-order MUSCL (Monotone Upstream-Centered Schemes for Conservation Laws) scheme was used for spatial interpolation. The PRESTO! (PREssure STaggering Option) scheme was selected as the discretization method for pressure. The PISO (Pressure-Implicit with Splitting of Operators) scheme was used for pressure–velocity coupling. Temporal discretization was a first-order implicit scheme. A dynamic mesh model built into FLUENT™ was employed to explicitly consider the sphere vibration.

To examine the effect that the fish body and fins have on the received dipole pressure signal, two body geometries were considered: body with fins extended (Figs 1 and 2) and body with fins retracted (this geometry is not shown). Also considered was a virtual-body case where the standard potential dipole source flow solutions (Pozrikidis, 1997) in an unbounded domain (without any internal fish boundaries) were evaluated at virtual lateral line locations identical to those on a real fish body. For all cases, the time step for integration was set at 1/100th of the vibration period. The 2-D flow field was initialized by simulating 200 time steps; the 3-D flow field was initialized by simulating 1000 time steps. This was done to allow transients to decay. After the initial start-up, flow velocity and pressure fields were saved at each time step for post-processing.

Fig. 1.

(A) 2-D computational domain for a circular cylinder (with cross-section identical to that of the sphere) vibrating nearby a sculpin cross-section. (B) A zoomed-in view of the meshes surrounding the cylinder and the sculpin. The vibrating motion of the cylinder was represented by a small rectangular deforming mesh zone that surrounds the immediate area of the cylinder. (C) Node locations around the sculpin cross-section where pressures at pore openings were evaluated. The convention for calculating the pressure difference is identified by the direction of the arrows. The pressure at the arrow tip (p1) is subtracted from the pressure at the arrow tail (p2). The distance between two consecutive pore openings is assumed to be 2 mm. dp/ds denotes the pressure gradient along the profile.

For the 2-D setup, pressure gradients were calculated along the sculpin cross-section by: Math(2) where p1 and p2 are pressure values at two consecutive nodal points on the profile, which were 2 mm apart by mesh construction and dp/ds denotes the pressure gradient along the profile. The ordering of p1 and p2 was defined by the arrows shown in Fig. 1C. For the 3-D setup, nodal points (Fig. 2C) were manually selected to approximately follow what might be considered the true locations of lateral line canals on the mottled sculpin [see fig. 1 in Coombs (Coombs, 2001)]. As the selected nodal locations were not evenly spaced, a spline curve for the pressure values was fitted to those locations for each of the lateral line sections (i.e. trunk canals, infraorbital canals, supraorbital canals, mandibular canals and occipital canal) (Fig. 2C). Next, pressure values were sampled every 2 mm on each fitted curve and then used in Eqn 2 to calculate pressure gradients. Also, for comparison purposes, the mid-plane cross-section of the 3-D fish body was identified, and pressure gradients were calculated along the mid-plane cross-section profile (Fig. 2D), similar to what was done for the 2-D setup.

Fig. 2.

(A) 3-D computational domain for a sphere vibrating nearby a sculpin body. (B) A zoomed-in view of the surface meshes of the sphere and of the sculpin, with nodes spaced approximately 1 mm apart on the sphere and 2 mm apart on the sculpin. The vibrating motion of the sphere was directly represented by a deforming mesh that surrounds the immediate vicinity of the sphere and is contained inside a small rectangular prism (not shown). (C) Node locations around the sculpin body where pressures were evaluated, including those points which fall directly on the sculpin's lateral line system (marked by various colors). (D) Points around the sculpin mid-plane cross-section, where pressures were interpolated from the pressure data computed on the 3-D sculpin surface. In both C and D, the convention for calculating the pressure difference is identified by the arrow direction. The pressure at the arrow tip (p1) is subtracted from the pressure at the arrow tail (p2). The distance between two consecutive pore openings is assumed to be 2 mm. dp/ds denotes the pressure gradient along the profile.

To validate the 3-D setup and simulation procedure, a few cases of a vibrating sphere above a flat plane wall were simulated. The mesh densities over the sphere surface and over the plane wall were the same as those used for the 3-D simulations that incorporated the fish body. The numerical results for the magnitudes of the pressure and pressure gradient (which are not affected by viscosity) showed excellent agreement with the PFT solution of a sphere vibrating above a plane wall (Fig. 3). The 2-D setup and simulation procedures were validated in a similar way (Rapo, 2009).

Flat-plate oscillatory boundary layer simulations and analytical model

3-D simulations were also performed of the formation of oscillatory boundary layers due to the vibration of the sphere. When a sphere vibrates near a stationary fish body in otherwise quiescent water, two oscillatory boundary layers will form, one on the sphere and one on the fish skin. Outside these two boundary layers, PFT is applicable (Fig. 3). The significance of the boundary layer on the sphere was recognized by Kalmijn (Kalmijn, 1988). A mathematical formulation was provided by van Netten (van Netten, 2006).

Both oscillatory boundary layers have the same Stokes viscous length scale,δ , defined as: Math(3) where ν is the kinematic viscosity of the fluid andω =2πf. (For f=50 Hz, which is a typical condition used in many experiments, δ∼80 μm.) Because of the oscillatory nature of the motion considered, viscous diffusion due to the presence of the wall cannot penetrate beyond a distance of order δ away from the wall [e.g. pp. 263–272 in Panton (Panton, 1996)]. The magnitude of the velocity of the oscillatory boundary layer of the fish skin is usually much weaker than that of the oscillatory boundary layer surrounding the sphere, with the difference in magnitude depending on the distance between the sphere and the fish's skin.

The 3-D computational mesh as described in the previous section was unable to resolve the oscillatory boundary layer along a fish-shaped body because the distances from centroids of the wall-adjacent control volumes to the fish skin were ∼290 μm, which is much larger than δ. To resolve this boundary layer, we would need to use a much finer near-wall mesh, which our current computational setup could not handle. To deal with this problem, we used CFD to calculate, instead, the flow along a flat plate produced by a nearby vibrating sphere. Anderson et al. (Anderson et al., 2001) have shown that flat-plate boundary layer theory can be used to model the boundary layer along a fish, and we adopted this assumption for the oscillatory case. The flat-plate geometry required fewer mesh points overall than a fish-shaped body while still allowing for a fine grid-point spacing at the wall (grid-point locations were at 0, 5, 13, 24, 39, 61 μm... away from the wall). The plane wall surface was covered by quadrilateral meshes. A deforming mesh zone, consisting of tetrahedral meshes, was present around the sphere so as to represent the motion explicitly. The volume between the deforming mesh zone and the boundary layer mesh zone close to the plane wall was divided into tetrahedral control volumes. The same numerical schemes and time step as previously described were used. Large parameter ranges as found in the literature were considered; velocity amplitude 7 mm s–1U0≤314 mm s–1, sphere radius 2.5 mm≤a≤18 mm and distance from sphere to fish body 1.1 cm≤r≤20 cm.

For the remainder of this section, we consider a sphere with vibration axis parallel to the wall. We define a Cartesian coordinate system such that the positive x-direction is parallel to the wall and to the axis of the sphere vibration, the positive y-direction is the wall-normal direction toward the sphere, the positive z-direction is chosen such that the defined xyz-coordinate system satisfies the right-hand rule and u, v and w are the x-, y- and z-velocity components, respectively. Using this Cartesian coordinate system, the strain rate tensor, S, can be defined as (e.g. Pozrikidis, 1997): Math(4) Also we define the overall strain rate (S) as: Math(5) S is a useful measure of the magnitude of fluid-parcel deformation irrespective of directional information. Within the oscillatory boundary layer at the fish's skin, bending of an SN cupula will be affected by the dynamic behavior of S surrounding the cupula as well as the whole velocity profile, u(y,t).

We use the CFD simulation of the oscillatory flat-plate boundary layer to validate a simpler analytical model, which we then use in the next section of this paper to analyze threshold velocity and S responses of SNs. The following is the analytical model. For a sphere vibrating parallel to a flat plane wall with r>>δ, the wall-parallel velocity profile along the wall-normal (y-) direction has an approximate solution: Math(6) This solution is for the flow created due to an oscillating plate, which is written in the frame of reference of the plate (e.g. Pozrikidis, 1997). The boundary-layer edge velocity is replaced by the wall slip flow velocity that corresponds to a potential dipole source placed at the center of the sphere and an image dipole source placed an equal distance below the wall to satisfy the potential flow `zero-normal-flow' wall boundary condition. C(a,δ) is the amplitude correction term due to the sphere boundary layer and ϕ is the phase delay. Both are calculated according to van Netten (van Netten, 2006) as: Math(7) Because this case involves only unidirectional flow, S is just| du/dy|, which is found by taking the y-derivative of the velocity profile (Eqn 6): Math(8) The velocity and strain rate profiles calculated using Eqns 6, 7, 8 are in excellent agreement with the results of CFD simulations for different parameter values (Fig. 4).

The maximum strain rate at the wall, Swall, is found by setting y=0 and taking the maximum, which occurs at (ωt+ϕ)=3π/4. This gives: Math(9) The maximum shear stress at the wall is just μSwall, where μ is the fluid dynamic viscosity, which is equal to 1.0×10–3 kg m–1 s–1. Eqn 9 is a useful measure of the forces acting to displace the cupula, especially when the height of the cupula, H, is unknown. If the height is known, another useful measure is maximum average strain rate, Saverage, along the cupula, which is defined as: Math(10) where Δu is the maximum velocity difference between the cupular tip and base over the whole vibration cycle.


Simulated prey-tracking sequence: the 2-D case

The initial three steps of a six-step, prey-tracking sequence recorded by Coombs and Conley (Coombs and Conley, 1997a) for a mottled sculpin were simulated using the 2-D setup (Fig. 5A) and the 3-D setup (Fig. 5B–D). The sphere had a radius of 3 mm, oscillating at a frequency of 50 Hz with source velocity amplitude of 0.18 m s–1 (peak-to-peak). These three positions were selected because they represented different relative locations between the fish's body (and pectoral fins) and the dipole source. In the starting location at the time of signal onset (Position No. 1, Fig. 5A,B), the sphere was slightly less than a body length away and was closer to the fish's tail than head. In the second position (Position No. 2, Fig. 5C), the sphere was also less than a body length away and was lateral to the point of pectoral fin insertion. In the third position (Position No. 3, Fig. 5D), the dipole source was in a more frontal location much closer to the head of the fish than the tail.

For 2-D CFD simulations, with the fish body present (columns 2 and 3 in Fig. 5A), the iso-pressure lines terminate on the fish body and are more concentrated around curved regions of the body and sharp edges of the fins. These local concentrations of pressure contours cannot be predicted from PFT without the fish body present (column 4 in Fig. 5A). When the fins are extended (column 2 in Fig. 5A), a zone of near constant pressure is created on the side of the fish closest to the source (the ipsilateral side) in a small, localized pocket behind the extended fin and side of the body. The presence of the fish body shields a large region of water from the dipole source, creating a region of near constant pressure along the contralateral side of the fish opposite the dipole source. Consequently, the pressure gradient seen by the lateral line on both the ipsi- and contralateral side tends toward zero as the fin insertion point (indicated by the red arrow in column 2 in Fig. 5A) is approached. This is in contrast to the case with a streamlined body present, where only the body curvature itself determines how the pressure field is affected (column 3 in Fig. 5A). In this case, there is no constant pressure region on the ipsilateral side of the fish near the pectoral fin. The iso-pressure lines of the hydrodynamic field without a body present are obviously undisturbed (column 4 in Fig. 5A).

Fig. 3.

Comparison between the 3-D computational fluid dynamics (CFD) simulations (symbols in C–F) and the potential flow theory (PFT) (solid lines in C–F) for a sphere vibrating above an infinite flat plane wall. The PFT solution consists of a dipole source above the wall and an image dipole source below the wall to satisfy the zero-normal-flow wall boundary condition. The axis of sphere vibration is either parallel (left column) or perpendicular (right column) to the wall. Iso-pressure contours of the instantaneous pressure field are depicted in (A) and (B). Line plots of the instantaneous pressure (C,D) and pressure gradient (E,F) along the wall are shown for three distances (1.8, 3.5 and 8.5 sphere diameters) between the sphere and the wall. U0, velocity amplitude of the sphere vibration; ρ, fluid density; ω=2πf; p1, p2, pressure values along the wall, which are used to calculate the pressure gradients, dp/ds.

Simulated prey-tracking sequence: the 3-D case

The same sculpin tracking sequence was simulated using the 3-D setup (Fig. 5B–D). The presence of the fish body perturbs the pressure field but to a lesser extent than for the 2-D case. Whereas there are clear differences between the results for fins extended and fins retracted in the 2-D case (columns 2 and 3 in Fig. 5A), there is no such detectable difference in the 3-D case (columns 2 and 3 in Fig. 5B). The reason is that, in the 3-D case, the dipole source flow field extends over and under the fins, as well as around them and therefore there is no longer any zone of near constant pressure behind the extended fins. In the 3-D case, the presence of the fish body still shields a region of water on the contralateral side from the dipole source but the overall effect of this shadow zone is weaker than for the 2-D case. The overall pressure field calculated from PFT without the fish body present is different from that obtained from the 3-D CFD simulations but the magnitude of the calculated pressure gradients along the virtual fish mid-plane profile are quite similar to those calculated from the 3-D CFD simulations. This is not the case in the 2-D simulation.

Fig. 4.

Comparison between the 3-D computational fluid dynamics (CFD) simulations (symbols) and the analytical solution (solid lines from Eqns 6, 7, 8) for the oscillatory boundary layer at the surface of an infinite flat plane wall, created by a sphere vibrating above and parallel to the wall. Vertical (y-) profiles are compared at the wall location directly below the equilibrium position of the sphere. The left column corresponds to the along-wall flow velocity, u, and the right column corresponds to the strain rate, S, and the shear rate, du/dy. A variety of sphere-to-wall distances (r), source vibration magnitudes (U0), frequencies (f), and sphere radii (a) are considered. Four examples are shown here: (A) f=30 Hz, r=10a U0=0.04 m s–1, phaseω t=0.4π, where ω=2πf and t the time. (B)f=45 Hz, r=20a, U0=0.1 m s–1, phase ωt=1.6π. (C)f=50 Hz, r=3.67a, U0=0.007 m s–1, phase ωt=2π. (D)f=75 Hz, r=40a, U0=0.03 m s–1, phase ωt=2π. δ, Stokes viscous length scale; U, fluid velocity just outside the oscillatory boundary layer.

The magnitude of the pressure gradient around the fish body is directly related to the location of the dipole source. Spatial variations are most striking when comparing the ipsilateral side of the fish closest with the dipole source and the contralateral side. However, equally noticeable are the differences in magnitude seen between the head and trunk of the fish. In the starting position (Fig. 5B), the magnitude of the pressure gradient is larger on the ipsilateral side near the end of the trunk (blue line) than at the front of the fish (green line). In the second position (Fig. 5C), the magnitude of the pressure gradient is similar for both the head and trunk (green line versus blue and red lines) but a clear difference still exists between the two sides of the fish (blue line versus red line). In the third position (Fig. 5D), the head of the fish is much closer to the dipole source, while the tail is further away. The pressure gradient amplitude has more than quadrupled for the front sections of the head canals (green line). Also, the contrast in pressure gradients between the two sides of the fish has widened (blue line versus red line).

Fig. 6 shows the pressure gradients along the lateral line canals of the mottled sculpin (Fig. 2C) based on the 3-D CFD results for the third position of the video sequence (Fig. 5D). Interestingly, the ipsilateral (positive numbers on the horizontal axis) pressure-gradient patterns along supraorbital, infraorbital and mandibular canals with different elevations (above and below the eye and along the lower jaw) but largely overlapping azimuths (rostro-caudal extents) tend to converge on the same, nearly redundant pattern. However, patterns on the ipsi- and contralateral side (negative numbers) are dramatically different. Ipsilateral patterns have steeper slopes and distinct zero-crossings (locations where the direction or sign of the pressure gradient changes from positive to negative). In this case, the true source location is near the zero-crossing point on the ipsilateral side.

Applying the oscillatory boundary layer solution to real experiments

We applied the analytical solutions Eqns 6 and 9 for the oscillatory boundary layer to the results of a number of previous neurophysiological and behavioral experiments designed to measure the threshold sensitivity of SNs in different species and under varying conditions. We used the reported experimental parameters to re-compute values of the velocity threshold at the tip of the SN cupula and the strain rate threshold at the wall, Swall. The parameters needed for the analysis are a, r, f, U0 and H. The following is a description of the different experiments used in the analysis and the results of the calculations.

Kroese et al. (Kroese et al., 1978) used a=1.55 mm oscillating at f=20 Hz parallel to the skin and perpendicular to the longitudinal axis of an SN of the clawed frog (Xenopus laevis), with U0=5×10–4 ms–1 and placed at a distance of r=3.75 mm from the skin. They assumed an SN cupula height of 100 μm and found a velocity threshold of 38μ ms–1 based on the potential flow equations. The Stokes viscous length scale for these particular experimental conditions calculated using Eqn 3 is δ=126 μm. Thus, the cupula would be fully immersed in the oscillatory boundary layer flow. From Eqn 6, we find that u at the tip of the cupula after being corrected for the boundary layer effects is 30 μms–1 (about 21% smaller than Kroese et al.'s threshold estimated from PFT). From Eqn 9, we find that the maximum Swall threshold is 0.45 s–1 and from Eqn 10 the maximum Saverage is 0.30 s–1.

Coombs and Janssen estimated velocity thresholds for SNs located along the trunk lateral line of mottled sculpin from neurophysiological measurements (Coombs and Janssen, 1990). The velocity field was produced by an oscillating sphere whose center was placed r=15 mm away from the fish trunk. The radius of the sphere was 3 mm and it was oscillated in a direction perpendicular to the substrate (up/down with respect to the fish) at frequencies in the range 10–500 Hz. The source velocity amplitude of the sphere corresponding to the threshold response, while not explicitly given, can be calculated using the dipole equations (e.g. Pozrikidis, 1997) and the results given in fig. 7 of Coombs and Janssen (Coombs and Janssen, 1990). We estimated that the peak-to-peak acceleration threshold (at the fish) for SNs at 10 Hz is –55 dB re. 1 m s–1 and that the threshold acceleration increases linearly by about 7.5 dB octave–1. From this, we determined that the corresponding source velocity amplitude of the sphere needed to produce the measured peak-to-peak acceleration in an unbounded fluid to be about 3.5 mm s–1 for 10 Hz, 5.3 mm s–1 for 50 Hz and 6.3 mm s–1 for 100 Hz. The corresponding velocity amplitude threshold at the tip of the SN cupula based on potential flow equations (with values doubled to account for the presence of the fish body) is 28 μms–1 at 10 Hz, 42μ ms–1 at 50 Hz and 50 μms–1 at 100 Hz. The velocity threshold values corrected for the oscillatory boundary layer effects (assuming H=100 μm) are 18 μms–1 at 10 Hz, 42 μms–1 at 50 Hz and 54 μms–1 at 100 Hz, where δ=176 μm at 10 Hz, δ=80 μm at 50 Hz andδ =56 μm at 100 Hz. The velocity threshold for the 10 Hz case is lower than the potential flow case due to the slowing of the flow near the wall, the velocity threshold for the 50 Hz case is unchanged and the velocity threshold for 100 Hz case is actually slightly higher than the potential flow case due to overshoot in the velocity profile (Fig. 7). {At certain distances from the wall the phase lag in viscous stresses is so great that the viscous and pressure terms actually add together. The combination of these forces accelerates the fluid to produce the overshoot [pp. 268–269 in Panton (Panton, 1996)].} The experimental parameter values give maximum Swall of 0.24 s–1 at 10 Hz, 0.78 s–1 at 50 Hz and 1.30 s–1 at 100 Hz. With H=100 μm, the maximum Saverage are 0.18 s–1 at 10 Hz, 0.42 s–1 at 50 Hz and 0.54 s–1 at 100 Hz.

Coombs et al. used a vibrating sphere to stimulate CNs and SNs and to see which subsystem was responsible for the orienting response in mottled sculpin (Coombs et al., 2001). When they disabled SNs, the fish was still able to orient towards the vibrating sphere. However, when they disabled CNs, response rates fell to spontaneous levels, indicating that the SNs were not used for orienting. An alternative explanation is that during the experiment, the stimulus to the SN was below its threshold velocity and strain rate. In their experiments, Coombs et al. used a sphere of radius of 3 mm placed 3–6 cm from the fish (Coombs et al., 2001). The source velocity was 5.5 cm s–1 for the 10 Hz signal and 9.0 cm s–1 for the 50 Hz signal. The maximum velocity at the tip of the cupula [taking viscous effects into account and assuming that H is 100 μm, as in Coombs and Janssen (Coombs and Janssen, 1990)] is 5–36 μms–1 for 10 Hz (range of values in all cases depend on the distance between the sphere and the fish) and is 11–89μ ms–1 for 50 Hz. The maximum Swall was 0.06–0.47 s–1 for 10 Hz and 0.21–1.66 s–1 for 50 Hz. With H=100 μm, the maximum Saverage are 0.05–0.36 s–1 for 10 Hz and 0.11–0.89 s–1 for 50 Hz. These values are above the SN threshold values measured by Coombs and Janssen (Coombs and Janssen, 1990) when the fish was 3 cm from the sphere but are below threshold values when the fish was 6 cm from the sphere.

Alternately, Abdel-Latif et al. showed that blind cave fish could orient toward a vibrating sphere even when the CN subsystem was disabled, indicating that the fish were relying on the SN subsystem to locate the sphere (Abdel-Latif et al., 1990). Their experimental parameters were a=2.5 mm, r=20 cm, f=10–90 Hz and displacement amplitudes = 0.2–1.4 mm. A positive behavioral response for frequencies 50 Hz and 70 Hz was reported for all amplitudes. The corresponding velocity thresholds (taking into account the oscillatory boundary layer effects) based on an SN height of 200 μm (Teyke, 1988) are 0.14–0.96 μms–1 at 50 Hz and 0.19–1.31μ ms–1 at 70 Hz (depending on the displacement amplitude). These values are well below the threshold values of the SNs of other species described above. In addition, the maximum Swall range is 0.0023–0.016 s–1 at 50 Hz and 0.0037–0.026 s–1 at 70 Hz, and the maximum Saverage range is 0.0007–0.048 s–1 at 50 Hz and 0.0010–0.0065 s–1 at 70 Hz. Again these are extremely low values, even acknowledging the possibility of signal startup transients, which from the present numerical simulation briefly increased the maximum wall strain rate up to 4 times the steady state value.


The results from the 2-D and 3-D simulations are significantly different. Qualitatively, the perturbation to the dipole pressure field due to the presence of fish body is much more prominent in the 2-D CFD results than in the 3-D CFD results. Similarly, the effect of the pectoral fin is also more visible in the 2-D than in the 3-D results. For our specific case, the 2-D and 3-D CFD simulations produced the `Mexican-hat' shaped pressure-gradient patterns with two zero-crossings for the head canal (green line of column 2 in Fig. 5D for 3-D, not shown for 2-D). However, the simulations produced completely different pressure-gradient patterns for the trunk canals (Fig. 5A versus 5B). Quantitatively, the most prominent difference is that the 2-D pressure-gradient magnitudes are 20 times larger than those calculated using the 3-D geometry. Taken together, this raises concerns about using 2-D results to calculate the hydrodynamic stimulus, at least for the present situation of a fish using its lateral line system to locate a vibrating sphere.

The 3-D results confirm that it is valid to use PFT to predict the spatial patterns of pressure gradient that affect the canal subsystem. The CFD simulations show that the fish body does perturb the dipole pressure field but only to a small extent, such that the pressure gradients evaluated along the lateral line locations are quite similar to those calculated from the 3-D virtual body case using PFT. This confirms the effectiveness of using PFT to estimate pressure-gradient patterns to the canal subsystem created by a vibrating sphere in still water (Goulet et al., 2008). The 3-D CFD results do show that the extended fins significantly distort the dipole pressure field locally. The fact that the lateral line trunk canals of the mottled sculpin are routed above the pectoral fins is probably of morphological significance, as it appears that this may limit the distortion in the received signal.

Fig. 5.

Comparison between 2-D calculations (A) and 3-D calculations (B), and comparison between 3-D computational fluid dynamics (CFD) (columns 2 and 3 in B–D) and 3-D potential flow theory (PFT, column 4 in B–D) results when the pectoral fins are either included (column 2) or excluded (columns 3 and 4). The calculations correspond to a video-recorded sequence of a mottled sculpin's step-by-step approach towards an artificial prey – in this case a sinusoidally vibrating sphere [from Coombs and Conley (Coombs and Conley, 1997a)]. The first three steps of the prey-tracking sequence are illustrated, including the initial orienting response at signal onset (A and B are for the same step) followed by two subsequent approach steps (C,D). Each block in columns 2–4 contains a plot of the iso-pressure contours of the normalized instantaneous pressure field that surrounds the sculpin and the vibrating sphere. Plotted in the lower panel are distributions of normalized pressure gradient along the three major portions of the lateral line canal system: the trunk canal on the ipsilateral (I) side of the body with respect to the dipole source (blue curve), the trunk canal on the contralateral (C) side of the body (red curve) and the frontal canals (F) on the head (green curve). T stands for the tail position. The pressure contours in column 4 correspond to the solution of a 3-D (2-D in A) dipole source in an unbounded fluid, unperturbed by the presence of the fish. The pressure gradient values plotted in column 4 are calculated from the unperturbed pressure field but the magnitudes are doubled to account for the potential flow wall-boundary condition of the fish and to facilitate comparison with the CFD solutions. The red arrow in column 2 of A indicates the pressure gradient results for the pectoral fin insertion points. U0, velocity amplitude; ρ, fluid density;ω =2πf; a, sphere radius; dp/ds, pressure gradient; p, pressure.

Fig. 6.

The complete dipole source pressure-gradient signal to the whole canal lateral line system, calculated from 3-D computational fluid dynamics (CFD) and evaluated for strike position No. 3 as shown in Fig. 5D. U0, velocity amplitude; ρ, fluid density;ω =2πf; dp/ds, pressure gradient.

Our higher resolved CFD boundary layer simulations using a flat-plate model highlight the fact that a vibrating sphere generates an oscillatory boundary layer at the fish's skin. The resulting flow velocity signal to the SNs and their responses are greatly affected by characteristics of the oscillatory boundary layer, which include its thickness (characterized by the Stokes viscous length scale) and the time-varying and height-varying flow magnitude and direction. Because of the presence of a boundary layer at the fish's skin, the velocity and velocity gradient patterns that act on the SNs cannot be predicted by PFT. In particular, there is a phase difference (Fig. 4A) between the flow velocity outside and inside the oscillatory boundary layer. Also, there is overshooting of the flow velocity inside the boundary layer (Fig. 4C,D), such that the velocity inside the boundary layer can be greater than the outside flow. Thus, the height of the cupula becomes an important parameter. For example, a cupula with a shorter height that is completely immersed in the boundary layer will experience a weaker integrated flow over its length than a taller cupula that extends outside the boundary layer. Consequently, a weaker response from the shorter cupula will be expected. Mukai et al. found that the feeding rate of the willow shiner larvae consuming Artemia was proportional to the cupular length on the larva body and saturated at lengths above ∼100 μm (Mukai et al., 1994). The trend of their fig. 1 curve relating the feeding rate to the cupular length corresponds well with the trend of the curves shown in Fig. 7 of the present study, which can be interpreted as the maximum tip-to-base velocity difference as a function of neuromast height for three different frequencies. Fig. 7 may provide an explanation to their results, provided that the larvae use the SNs to detect appendage-beating movement of Atemia.

There is a minimum cupular displacement that must occur in order for the neuromast to respond, corresponding to a minimum velocity whose drag force causes the displacement. This is called the velocity threshold. Although the velocity threshold may be species-specific, the CNs probably respond to internal canal fluid velocities as low as 1–10 μms–1 (van Netten, 2006). Inside the subdermal canal, the velocity is driven by the pressure difference between pore openings, which translates to an acceleration threshold of 0.1–1 mm s–2 (van Netten, 2006). The SNs respond to as little as 25–60μ ms–1 (for a review, see van Netten, 2006). These results were obtained using PFT with the assumption that the cupular height is larger than the boundary layer thickness. Using our analytical treatment, we showed that the oscillatory boundary layer formed at the fish's skin should be taken into account when estimating these velocity thresholds. The oscillatory boundary layer correction can be significant depending on the cupular height relative to the boundary layer thickness. Re-analysis of Coombs and Janssen (Coombs and Janssen, 1990) gives a threshold velocity range of 18–54 μms–1 over a frequency range of 10–100 Hz, when the oscillatory boundary layer effects are taken into account. The velocity threshold can also be defined as the maximum tip-to-base velocity difference experienced by the cupula, because the base velocity is zero. This is species-specific and depends on a number of parameters, including cupular height and shape, as well as internal cupula stiffness (McHenry et al., 2008).

From CFD outputs, one may calculate maximum wall-strain rates, which are a measure of the near-wall fluid-parcel deformation probably experienced by the SNs and are independent of cupular height. For the mottled sculpin of the Coombs et al. (Coombs et al., 2001) experiment, the maximum wall strain rates given in the previous section corresponded to the threshold values of SNs (either just above or just below, depending on the starting distance of the fish from the vibrating dipole source). At the same time, the acceleration experienced by the lateral line canal subsystem (based on Fig. 6) is 5–50 mm s–2, which is many times stronger than the reported acceleration detection thresholds of CNs. This is consistent with the conclusion that the acceleration-responsive CNs, rather than the velocity-responsive SNs, mediate the mottled sculpin's orienting behavior to the vibrating sphere (e.g. Coombs et al., 2001). However, all the SN thresholds calculated for the blind cave fish of the Abdel-Latif et al. (Abdel-Latif et al., 1990) experiment are extremely low as compared with those numbers obtained for other experiments. This raises concerns regarding claims that the SNs were used for detecting the dipole source, especially since these otophysan fish have pressure-sensitive ears that could have easily detected the pressure changes of a nearby dipole source of the same frequency (Montgomery et al., 2001).

Fig. 7.

The maximum tip velocity experienced by the superficial neuromasts (SNs) as a function of the neuromast height. The plot was generated analytically based on the experimental conditions used by Coombs and Janssen (Coombs and Janssen, 1990)– sphere radius, r, of 3 mm vibrating in a direction perpendicular to the frontal plane of the fish at frequencies, f, of 10, 50 and 100 Hz and with velocity amplitude, U0, of 3.5 mm s–1, 5.3 mm s–1 and 6.3 mm s–1. The distance from the sphere center to the fish surface was 15 mm. The solid circle is the potential flow result for an SN with cupula height of 100 μm. u, x velocity component.


sphere radius
C, C1, C2
amplitude correction term
computational fluid dynamics
canal neuromast
pressure gradient
shear rate
mathematical constant e
vibration frequency
cupula height
p1, p2
pressure values at two consecutive nodal points on the profile
potential flow theory
distance from the sphere to the fish body
strain rate tensor
strain rate
average strain rate
superficial neuromast
wall strain rate
x velocity component
velocity profile
velocity amplitude of the sphere vibration
fluid velocity just outside the oscillatory boundary layer
y velocity component
z velocity component
Stokes viscous length scale
maximum velocity difference between cupular tip and base
fluid dynamic viscosity
fluid kinematic viscosity
fluid density
phase delay


The authors are grateful to the Woods Hole Oceanographic Institution Academic Programs Office for graduate research support provided to M.A.R. The authors gratefully acknowledge NSF for support of this research under grants IOS-0718506 to H.J. and IBN-0114148 to M.A.G. The authors are grateful to two anonymous reviewers, who provided insightful comments that improved the final paper. The authors thank A. H. Techet for helpful discussions.


View Abstract