## Abstract

In Part I of this two-part study, the coupled flows external and internal to the fish lateral line trunk canal were consecutively calculated by solving the Navier–Stokes (N–S) equations numerically in each domain. With the external flow known, the solution for the internal flow was obtained using a parallelepiped to simulate the neuromast cupula present between a pair of consecutive pores, allowing the calculation of the drag force acting on the neuromast cupula. While physically rigorous and accurate, the numerical approach is tedious and inefficient since it does not readily reveal the parameter dependencies of the drag force. In Part II of this work we present an analytically based physical–mathematical model for rapidly calculating the drag force acting on a neuromast cupula. The cupula is well approximated as an immobile sphere located inside a tube-shaped canal segment of circular cross section containing a constant property fluid in a steady-periodic oscillating state of motion. The analytical expression derived for the dimensionless drag force is of the form , where |*F*_{N}| is the amplitude of the drag force; |*P*_{L}−*P*_{R}| is the amplitude of the pressure difference driving the flow in the interpore tube segment; *d*/*D* is the ratio of sphere diameter to tube diameter; *L*_{t}/*D* is the ratio of interpore tube segment length to tube diameter; and is the oscillating flow kinetic Reynolds number (a dimensionless frequency). Present results show that the dimensionless drag force amplitude increases with decreasing *L*_{t}/*D* and maximizes in the range 0.65≤*d*/*D*≤0.85, depending on the values of *L*_{t}/*D* and . It is also found that in the biologically relevant range of dimensionless frequencies 1≤≤20 and segment lengths 4≤*L*_{t}/*D*≤16, the sphere tube (neuromast–canal) system acts as a low-pass filter for values *d*/*D*≤0.75, approximately. For larger values of *d*/*D* the system is equally sensitive to all frequencies, but the drag force is significantly decreased. Comparisons with N–S calculations of the drag force show good agreement with the analytical model results. By revealing the parameter dependencies of the drag force, the model serves to guide biological understanding and the optimized design of corresponding bioinspired artificial sensors.

## 1. Introduction

In Part I of this two-part study, the nature, function and hydrodynamic mode of the operation of the lateral line trunk canal (LLTC) were described by reference to prior work and results from a numerical model. Briefly, the LLTC consists of a long subepidermal tube-like conduit aligned along the length of the sides of a fish. It has pore-like openings periodically distributed along its length which connect the fluid external to the canal to the fluid internal to it. Between each pair of pores, inside the canal, there is a neuromast attached to the inner tube wall. Fluid motions external to the LLTC create spatial–temporal distributions of pressure, which are impressed at the pore openings and propagated into the LLTC. In this way, the external time-varying pressure gradients arising between pairs of pores produce corresponding time-varying flows (between pairs of pores) inside the LLTC. In turn, these flows produce drag forces that act on the neuromasts1 (between pairs of pores), thus causing nanometre-scale displacements of them. The displacements activate hair cells inside the neuromasts that inform the fish concerning the magnitude and direction of the flow internal to the canal. The information obtained from many neuromasts distributed along the length of the LLTC allows the fish to determine important characteristics of the external surrounding flow.

The results presented in Barbier & Humphrey (2009) are based on numerical solutions of the Navier–Stokes (N–S) equations for the flows external and internal to the LLTC. In their approach, the unsteady flow field external to the LLTC is solved for first. In the case investigated, this consists of the vortical flow in the wake of a prism-like object moving at the same speed as, but upstream and to one side of, a pikeperch (*Sander lucioperca*). In this one-way coupled system, the pressure field associated with the external flow passing over the pikeperch drives the flow inside its LLTC. This internal flow is then solved for between a consecutive pair of pores with a neuromast present halfway between them, the solution yielding the drag force acting on the neuromast as a function of the relevant problem parameters. While this approach can be made to be as physically exact and numerically accurate as desired, this comes at significant computational costs. This is for three reasons: (i) to accurately resolve the internal flow it is necessary to use refined grids capable of capturing the complex neuromast–canal geometry, (ii) each calculation case performed corresponds to one combination of the relevant dimensionless parameters and requires large computer time and storage, and (iii) the resolution of parameter dependency trends requires a large number of individual calculation cases. Such a problem solving approach is tedious and inefficient. As a consequence, and owing to the relative simplicity of the flow internal to the LLTC relative to the flow external to it, we ask: given the spatial–temporal characteristics of the flow external to the LLTC, in particular, the pressure field, is there a generally applicable, analytically based, physical–mathematical model capable of quickly and accurately predicting the parameter dependencies of the drag force acting on a neuromast?

The goal of this study has been to develop such a predictive theoretical model; one capable of capturing the *essential physics* determining the drag force induced on a neuromast inside the LLTC as a consequence of the external–internal flow coupling described above. As a by-product, the model should yield the purely physical dimensionless parameter dependencies of the drag force. With such a model in hand, it is then possible to explore the following two interrelated questions.

How does the drag force acting on a neuromast inside the LLTC vary as a function of the parameters affecting it? In particular, is there a ‘purely physical’ optimal neuromast size?

What are some of the possible implications to biology, and to the design of corresponding bioinspired neuromast-like flow motion sensors, of the drag force parameter dependencies observed?

## 2. The physical model

### 2.1 Geometry and pressure conditions

With reference to figure 1, we consider a segment of the LLTC bounded by two pores. The segment has the shape of a tube of length *L*_{t}=2*L*+*d* and constant diameter *D*. Halfway along the tube is a neuromast shaped as half a prolate spheroid of base diameter *d*_{N} and height *h*_{N}. For estimating the drag force acting on the spheroidally shaped neuromast, it is especially convenient to approximate the neuromast as a sphere of geometrical mean diameter . The goodness of this approximation is discussed in detail in §4. The approximation has the virtue of yielding the volume of the half-prolate spheroid exactly and its surface area to within 0–20 per cent depending on the value of *h*_{N}/*d*_{N} in the range 0.75≤*h*_{N}/*d*_{N}≤4. Thus, the spherical shape is expected to yield good estimates of the viscous and added mass contributions to the drag force acting on the neuromast.

Again with reference to figure 1, the pressure of the external flow impresses itself upon the left and right pore openings straddling the neuromast. Owing to the relatively short lengths of the pore openings, it is assumed that this is the same pressure difference that drives the flow inside the segment. The pressure drop inside the segment is due to the wall friction in the two sub-segments of length *L*, and to the wall friction plus the drag caused by the neuromast in the sub-segment of length *d*. Thus, we take(2.1)

Because, as shown below, the Reynolds number of the flow in the segment is much smaller than unity it follows that (*P*_{L}−*P*_{l})=(*P*_{r}−*P*_{R}) and, therefore,(2.2)In the biological situation, the quantity *P*_{L}−*P*_{R} generally consists of a constant value with an oscillatory component superimposed. The constant value is trivial to deal with and relatively uninteresting. Therefore, we focus attention on the oscillatory component and assume a cosinusoidal pressure difference waveform given by(2.3)where |*P*_{L}−*P*_{R}| and *ω* (≡2*πf*) are the wave amplitude and angular frequency, respectively. The oscillating pressure difference *P*_{L}−*P*_{R} sets up a corresponding oscillating flow inside each of the three tube sub-segments.

### 2.2 Pressure-driven oscillating laminar flow in a tube

The pressure-driven oscillating laminar flow inside a tube of diameter *D* has been solved analytically and results are reported in Rosenhead (1963), Goldstein (1965), White (1974) and Schlichting (1979). The dimensionless parameter governing such a flow is given by the kinetic Reynolds number , where *ν* (≡*μ*/*ρ*) is the kinematic viscosity of the fluid, *μ* being the dynamic viscosity and *ρ* the density. The final solution for the oscillating fluid velocity in the tube (for values <2000, approx., to ensure laminar flow) depends in a complicated way on the Bessel function of the first kind with imaginary arguments. However, the solution admits simpler forms for very small and very large values of . For very small values of , the bulk of the oscillating flow in the tube is governed primarily by a balance between pressure gradient and viscous forces, whereas for very large values of , it is governed primarily by a balance between pressure gradient and inertial forces. Characteristic velocities, *V*_{ch}, for the fluid near the tube centre-line are given by White (1974).

Viscous-dominated flow (≪4)(2.4)In this case, the characteristic velocity is in phase with the pressure oscillation.

Inertia-dominated flow (≫4)(2.5)In this case, the characteristic velocity lags the pressure oscillation by 90 degrees.

In equations (2.4) and (2.5) the quantity *K* is related to the pressure gradient, d*P*/d*x*, inside the tube by(2.6)Thus, for a tube of length *L*_{t} with an imposed pressure difference given by equation (2.3) it follows that *K*=|*P*_{L}−*P*_{R}|/(*ρL*_{t}).

In what follows we relax the above conditions to require simply that <4 for equation (2.4) to apply and >4 for equation (2.5) to apply. For the analytical model this implies that viscous ↔ inertia transitions for the flow in the tube occur at ≈4 when, in fact, they take place gradually over a range of centred about ≈4. As will be shown, in the presence of a spherically shaped neuromast in the tube segment the approximation is good for values of *d*/*D*≤0.75 and yields correct qualitative trends for *d*/*D*>0.75.

### 2.3 Evaluation of the pressure drops in the three tube sub-segments

Equations (2.4)–(2.6) can be used to obtain expressions for (*P*_{L}−*P*_{l}), (*P*_{l}−*P*_{r}) and (*P*_{r}−*P*_{R}) in terms of the characteristic velocities *V*_{L} and *V*_{d} in the respective tube sub-segments (figure 1).

Using equation (2.6) in the two tube sub-segments of length *L*, we have(2.7)and, using equations (2.4) and (2.5), we get the following.

Viscous-dominated flow in sub-segments *L* (<4)(2.8)

Inertia-dominated flow in sub-segments *L* (>4)(2.9)

Using equation (2.6) in the tube sub-segment of length *d*, we have(2.10)and, using equations (2.4) and (2.5), we get the following.

Viscous-dominated flow in sub-segment *d* (<4)(2.11)

Inertia-dominated flow in sub-segment *d* (>4)(2.12)

In these last two cases, the kinetic Reynolds number is given by where the quantity *D*_{H} (≡*D*−*d*) is the hydraulic diameter (White 1974) of the cross section of the tube defined by the space between the spherically approximated neuromast and the tube. From mass conservation of the flow between tube sub-segments, it follows that(2.13)

### 2.4 Evaluation of the characteristic velocity in the tube sub-segment of length *d*

<4 in sub-segments *L* and <4 in sub-segment *d*.

Combining equations (2.2), (2.8), (2.11), (2.13) and the definition of *D*_{H} and solving for *V*_{d} yields(2.14)

>4 in sub-segments *L* and >4 in sub-segment *d*.

Combining equations (2.2), (2.9), (2.12), (2.13) and the definition of *D*_{H} and solving for *V*_{d} yields(2.15)Equation (2.14) applies to flows for which *both* and are *less than 4* in the respective tube sub-segments (viscous-dominated flows). Equation (2.15) applies to flows for which *both* and are *greater than 4* in the respective tube sub-segments (inertia-dominated flows). Because *D*_{H}≤*D*, it is clear that if <4 then <4 also and, for the same reason, if >4 then >4 also. While it is not possible for <4 and >4, it is possible for >4 and <4. This case corresponds to a flow that is inertia-dominated in the tube sub-segments of length *L* and viscous-dominated in the tube sub-segment of length *d*. We examine this final case next.

>4 in sub-segment *L* and <4 in sub-segment *d*.

Combining equations (2.2), (2.9), (2.11), (2.13) and the definition of *D*_{H} and solving for *V*_{d} yields(2.16)In obtaining equation (2.16) it is tacitly assumed that, to a good approximation, the time variation of the characteristic velocity in sub-segment *d* is dictated by that in sub-segment *L*. Equations (2.14) and (2.16) predict realistic values of *V*_{d} as a function of *d*/*D*. By contrast, equation (2.15) yields finite values of *V*_{d} for *d*/*D*=1, which is physically unrealistic. We show below that: (i) there exists a critical value of below which equation (2.15) applies and above which equation (2.16) should be used instead, and (ii) this critical value of is such that, in practice, equation (2.16) applies for all cases with >4 of interest to this study.

## 3. Evaluation of the drag force acting on a neuromast inside the LLTC

### 3.1 The drag force equation

An equation for the drag force, , acting on a sphere of diameter *d* oscillating cosinusoidally at frequency *ω* in a still fluid of density *ρ* and dynamic viscosity *μ* has been derived by, among others, Lamb (1932) and Landau & Lifshitz (1986). The derivation of this equation is premised on the assumption that the nonlinear convective (or spatial) acceleration terms in the momentum equation for the fluid are negligible. For this to be true at low oscillation frequencies, for which , it is necessary that where *a* is the amplitude of the oscillation. In turn, for this to be true at high oscillation frequencies for which , it is necessary that *a*/*d*≪1. Calling *V*_{S} the oscillating sphere velocity and defining the characteristic Reynolds number of the sphere by *Re*_{d} (≡*d*|*V*_{S}|/*ν*, with |*V*_{S}| the amplitude of *V*_{S}), and noting that *aω*≈|*V*_{S}|, the condition *ωad/ν*≪1 at low frequencies is equivalent to *Re*_{d}≪1. For the cases of interest here this condition is satisfied at all frequencies for the flow past a neuromast inside the LLTC.2 The additional condition that *a*/*d*≪1 at high frequencies is equivalent to stating |*V*_{S}|/(*ωd*)≪1. This condition can be rewritten as and is also satisfied for the ranges of *d*/*D* and of interest to this study. Note that, in contrast to low frequencies, at high frequencies *Re*_{d} need not be small.

Subject to the above constraints, the drag force is given by(3.1)Equation (3.1) is written in its most fundamental form such that the quantities *δ* and *d* appear naturally together as a dimensionless ratio. The first term on the right-hand side represents the dissipative component of the drag force acting on the sphere due to form and friction drag. The second term corresponds to the inertial component of the drag force in potential flow, associated with a volume of fluid, or ‘added mass’, which must be accelerated with the sphere. In these two terms, the quantity *δ* is a characteristic viscous length scale given by(3.2)

Using theoretical analyses provided by, for example, Batchelor (1970) and Brennen (1982), equation (3.1) can be converted to yield the drag force acting on a fixed sphere in an oscillating fluid. This is accomplished by noting that there is an additional contribution to the total force acting on the sphere, referred to by Batchelor (1970) as an effective ‘buoyancy’ force, ‘analogous to the force on the body arising from the action of gravity on the fluid’. The result is(3.3)where(3.4)is the far-field oscillating fluid velocity written in cosinusoidal form. This approach for obtaining the drag force acting on a spherically shaped neuromast simplifies the problem considerably while retaining all its essential physics (see van Netten 2006).

Equation (3.3) yields the drag force acting on a fixed sphere in a flow oscillating freely far from the sphere. By contrast, depending on the frequency, the drag force acting on a fixed sphere due to an oscillating fluid inside a tube may be affected by the tube wall. When such tube wall–sphere interactions arise, the expectation is that the drag force acting on the sphere will be larger than for the freely oscillating flow condition (Brennen 1982). For the ranges of *d*/*D* and of interest to this study, the flow in the tube sub-segment containing the spherically approximated neuromast is generally viscous dominated (<4). Therefore, tube wall–sphere interactions will be small if the thickness of the boundary layer around the sphere is less than the distance between the sphere surface and the tube wall, that is if . Using equation (3.2) allows this constraint to be written as . Thus, tube wall–sphere interactions will be small and equation (3.3) will be applicable if the kinetic Reynolds number satisfies the latter constraint. For the conditions of this study, this is the case for *d*/*D*<0.1, approximately.

However, it is possible to modify equation (3.3) to apply to the larger values of *d*/*D* by recognizing that when wall tube–sphere interactions arise the velocity gradient at the sphere surface is larger relative to the freely oscillating flow condition. The increased velocity gradient is due to the increased speed of the flow passing over the sphere obstructing the tube which, effectively, works to reduce the thickness of the boundary layer around the sphere. In this case, a more appropriate measure of the boundary-layer thickness is given by which ensures that the boundary layer around the sphere does not extend beyond the tube wall. This modification retains the essential problem physics and, as shown below, yields results in good agreement with calculations based on numerical solutions of the N–S equations.

The *practical conclusion* is that equation (3.3) can be used by setting when , and by setting when . In equation (3.13) below, this is equivalent to setting when , and when , in the corresponding *d*/*δ* terms.3

### 3.2 The drag force amplitude and phase angle

Calling(3.5a)and(3.5b)equation (3.3) reads(3.6)or, equivalently, using equation (3.4)4 for *V*_{F},(3.7)Equation (3.7) can be written in the form(3.8)where(3.9a)is the amplitude of *F*_{S} and *β* is the phase angle given by(3.9b)We now use equations (3.5*a*) and (3.5*b*) and equation (3.9*a*) to obtain the amplitude of the drag force |*F*_{N}| acting on a spherically shaped neuromast due to an oscillating flow inside a tube-shaped segment of a LLTC. Of interest are cases 2.7, 2.8 and 2.9 above for which equations (2.14)–(2.16), respectively, provide the characteristic velocities corresponding to equation (3.4) above. Thus, for each case, the quantity *C* in equation (3.9*a*) is known. After algebraic manipulation the final expressions for |*F*_{N}| in dimensionless form are given by the following.

<4 in sub-segments *L* and <4 in sub-segment *d*(3.10)

>4 in sub-segments *L* and >4 in sub-segment *d*(3.11)

>4 in sub-segments *L* and <4 in sub-segment *d*(3.12)

In these equations, the quantity is given by5(3.13)

By reference to figure 1 we see that *L*=(*L*_{t}−*d*)/2 and it follows that equations (3.10)–(3.12) are all of the functional form(3.14a)thus revealing the dimensionless parameter dependencies of the (dimensionless) drag force acting on a spherically shaped neuromast in a tube-shaped LLTC. Note that equations (3.10) and (3.12) yield the same analytical solution for the special case when =4.Algebraic manipulation of equation (3.9*b*) allows the drag force phase angle to be written as(3.14b)This expression shows that *β* is a function of *d*/*D* and only, the dimensionless tube length *L*_{t}/*D* playing no role. For the same reasons explained above in relation to equation (3.13), in equation (3.14*b*) the evaluation of in the internal parenthetical term in the numerator, and in the term in the denominator, requires setting when , and when .

## 4. Constraints affecting the analytical model and its accuracy

Before generating and interpreting calculations based on the analysis presented above, it is appropriate to examine the constraints embedded in the modelling approach due to the physical approximations made and the expected influence of these limitations on the accuracy of the analytical results. Two types of constraints affect the analysis depending on whether they apply to the flow in the tube approximating the canal segment or to the flow over the sphere approximating the neuromast. These are discussed in turn. The section concludes with an evaluation of the accuracy of the analytical model by reference to numerical solutions of the problem of interest based on the N–S equations.

### 4.1 Constraints affecting the flow in the tube

The analysis of the flow in the tube approximating a segment of the LLTC assumes developed, steady-periodic motion at low Reynolds number (

*Re*_{D}≡*D*|*V*_{L}|/*ν*).6 As in Part I, here we also find that this Reynolds number is*Re*_{D}≤10^{−3}and the laminar nature of the flow is not in question. However, the analysis is mute concerning initial transients in the motion starting from rest. This is not a serious limitation as the temporal development of the flow is expected to occur in a few cycles. Numerical solutions based on the N–S equations for the flow past a simulated neuromast in a segment of a LLTC, presented in Part I as well as in the present study, support this. In addition, the following scaling argument applies. Since the flow in the bulk of the tube is established by viscous diffusion of momentum in the radial direction, a typical time scale to build up the wall boundary layer is given by . Calling*T*=2*π*/*ω*the period of flow oscillation, we require that*τ*/*T*<1. Making substitutions yields . In this study 0.16≤≤3.18, and we conclude that the temporal development of the flow is established within less than 3.2 cycles, approximately.To ensure that the flow in the tube is not adversely affected by fluid turning 90 degrees to enter and leave the tube at the pore locations, it is necessary to show that the flow at these locations develops spatially over a very short length,

*L*_{dev}. For an estimate of this length we consider steady laminar flow in a tube, for which the dimensionless hydrodynamic development length is given by*L*_{dev}/*D*=0.06*Re*_{D}(Fox*et al*. 2008). Since*Re*_{D}≤10^{−3}it follows that*L*_{dev}/*D*≤6×10^{−5}, approximately, and the condition is readily satisfied. Additional support for the shortness of the development length is provided in fig. 7 of Part I.It has been assumed that the flow passing through the three tube sub-segments (figure 1) incurs the pressure drop given by equation (2.1), and that each of the terms in equation (2.1) can be separately evaluated from the analysis for the unobstructed oscillating flow in a tube. This implies that any additional pressure drops associated with spatial accelerations and decelerations of the flow entering and leaving the narrower central tube sub-segment are negligible. The data in fig. 16 of Part I support this, but for a rigorous justification it is necessary to show that spatial accelerations of the flow are small relative to viscous terms at low values of the kinetic Reynolds number, meaning for <4, and small also relative to temporal acceleration terms at high values of the kinetic Reynolds number, meaning for >4. Using the scaled approximations , and together with equation (2.13) where, again,

*T*=2*π*/*ω*, yields the constraints: for <4; and for >4. Both constraints are readily satisfied in the ranges 0.01≤*d*/*D*≤0.99 and 1≤≤20 of interest to this study.

The analysis assumes that the transition between viscous-dominated and inertia-dominated flow regimes in a tube-like canal segment containing a neuromast occurs at a kinetic Reynolds number ≈4. We know from theoretical results for unobstructed tube flows, and from our own N–S based numerical calculations with a neuromast present in the tube, that such a transition occurs more gradually over a range of near to 4. Notwithstanding, the approximation is expected to work well for tube flows with relatively small neuromasts, and comparisons between results from the analytical model and numerical solutions based on the N–S equations show that it is good for spherically shaped neuromasts with

*d*/*D*≤0.75 and yields qualitatively useful results for the larger values of*d*/*D*.Based on the N–S numerical results of Part I, the present analysis neglects the pressure drop associated with the length of a pore connecting the external and internal flows. It is assumed that the pore length is short enough that the external pressure imposes itself through the pore directly on the flow inside the canal. In those cases where this is not the case, the effect of pore length can be included in the analysis, but may require considering the effects of adjacent pores by solving a tube network. Also as in Part I of this study, it is assumed that the pore diameter and canal diameter scale within a factor ranging between 0.2 and 2. While investigations of pore length and diameter effects have been beyond the scope of this study, none of the parameter dependencies predicted here is significantly altered by either of these two factors.

The specification of the right-hand side of equation (3.14

*a*) yields the left. However, to calculate a*dimensional*value for the drag force amplitude |*F*_{N}| requires that*D*and |*P*_{L}−*P*_{R}| be known. In the biological situation, the amplitude of the pressure difference, |*P*_{L}−*P*_{R}|, is the result of flow oscillations external to the body surface of the fish and must be evaluated experimentally, or computationally as in Part I.

### 4.2 Constraints affecting the flow over the sphere

It is assumed that the spherically approximated neuromast is immobile relative to the fluid moving past it in the tube. In reality a LLTC neuromast is displaced by a few tens of nanometres whereas the fluid around it is displaced by a few tens of micrometres. Thus, the assumption of an immobile sphere is acceptable for purposes of evaluating the drag force acting on it.

Prior to making the spherical assumption, the LLTC neuromast was approximated as having the shape of half a prolate spheroid of height

*h*_{N}and base diameter*d*_{N}(figure 1). An even more general approximation would allow the base of the spheroid to be elliptical in shape rather than circular. However, because we do not expect the ratio of the large to small axes of such an elliptic base to be much larger than 2 for a typical neuromast, attention has been restricted to the mathematically convenient half-prolate spheroid approximated as a sphere. It is readily shown that a sphere with geometrical mean diameter given by has a volume exactly that of the half-prolate spheroid and a surface area that matches that of the half-prolate spheroid to within 1.5–20 per cent depending on the value of*h*_{N}/*d*_{N}in the range 0.75≤*h*_{N}/*d*_{N}≤4. The volume and surface area of the half-prolate spheroid are given by (Beyer 1986)(4.1a)(4.1b)with the eccentricity,*ϵ*, equal to(4.2)The corresponding volume and surface area of a sphere with geometrical mean diameter are given by(4.3a)(4.3b)From these results it is clear that a sphere with diameter yields the volume of a half-prolate spheroid exactly. The ratio

*S*_{S}/*S*_{HPS}is plotted as a function of*h*_{N}/*d*_{N}in figure 2 and the result shows that the surface of such a sphere approximates that of the half-prolate spheroid to within an error less than ±20 per cent (and as small as 0%) in the range 0.75≤*h*_{N}/*d*_{N}<4. From this analysis we may also safely conclude that the volume and surface area of a more generally shaped neuromast, described by a half-prolate spheroid of height*h*_{N}, length*l*_{N}and width*w*_{N}(∼*l*_{N}), are well approximated by those of a sphere with geometrical mean diameter .

### 4.3 Choice of model equation for calculating the drag force

In performing a model-based calculation of the drag force acting on a neuromast, it is necessary to establish which of equations (3.10), (3.11) or (3.12) applies, and over what range. Because the ratio is given by(4.4)fixing any two dimensionless quantities in this equation fixes the third. Figure 3 shows a plot of as a function of *d*/*D*. If, for example, a value of *d*/*D*=0.375 is fixed, then =0.39 and, among many, the following pairs of (, ) values are possible: (0.39, 1), (6.24, 16) and (3.12, 8). In the first case (case 3.1) both and are less than 4 and equation (3.10) applies. In the second case (case 3.2) both and are larger than 4 and equation (3.11) applies. In the third case (case 3.3) >4 and <4 and equation (3.12) applies. In practice, both cases 3.2 and 3.3 are accurately predicted by equation (3.12) over the whole range of *d*/*D*. That this is the case is illustrated by the following argument. We note from equation (4.4) or figure 3 that for every there is a value of *d*/*D* below which >4 and above which <4. In the latter case equation (3.12) must be used to calculate the drag force, and in the former equation (3.11), in principle, should be used. The *d*/*D* values for which passes through the value 4 are given by *d*/*D*=0, 0.293, 0.423, 0.500 for =4, 8, 12 and 16, respectively. For each of these the corresponding critical value of *d*/*D* fixes the point above which equation (3.11) ceases to be applicable (the drag force increasing unrealistically with increasing *d*/*D*) and equation (3.12) must be used; see also figure 9 discussed below. The outcome is that, for all practical purposes, equation (3.12) can be used for all oscillating flows with >4 (cases 3.2 and 3.3) over the whole range of *d*/*D*, with very small error for small values of *d*/*D* where, in principle, equation (3.11) is more accurate.

### 4.4 Assessment of the accuracy of the analytical model

The physical accuracy of the analytical model has been assessed by reference to N–S based calculations of the drag force acting on a neuromast due to an oscillating flow in a tube-like canal segment. Such solutions were presented in Part I; see fig. 11 of that paper for the case of a neuromast shaped as a parallelepiped. Additional solutions are provided here for oscillating flows in tube segments containing spherically approximated neuromasts. For this, full three-dimensional, unsteady forms of the N–S equations are solved for the primitive variables velocity and pressure. The no-slip velocity condition is specified at all solid surfaces. The flow field in the tube segment is driven by a sinusoidally time varying pressure difference. To this end, the pressure at one pore is kept at a constant (zero) reference value while the pressure at the other is varied according to equation (2.3). The drag force acting on a spherically shaped neuromast is determined by surface integration of the combined pressure and shear stress contributions to it. For both the analytical and numerical calculations of this study, the fluid in the tube is taken to be water with physical properties fixed at 21°C (provided below).

The additional N–S based numerical calculations presented here were obtained with CFX, an algebraic, multigrid element based, finite volume commercial code. CFX uses a second-order Euler method with inner iterations set to achieve a preset tolerance criterion within each time step, or with a maximum number of iterations set per time step. The procedure employs a full hexahedral mesh using ICEM-CFD, a commercial meshing tool. An O-grid is used to define the field around the sphere in the tube in order to prescribe a fully hexahedral mesh.

In the numerical calculations of this study, the diameter of the tube segment is set to 0.25 mm and two lengths (2 and 4 mm) are examined yielding *L*_{t}/*D*=8 and 16. For each of these cases three sphere diameters are set (0.0625, 0.125 and 0.1875 mm) corresponding to *d*/*D*=0.25, 0.50 and 0.75. For each of these cases calculations are performed at flow oscillation frequencies of 10, 25, 40, 60, 80, 100, 150 and 200 Hz corresponding to kinetic Reynolds numbers =1, 2.5, 4, 6, 8, 10, 15 and 20. All calculations are started from a flow initially at rest and taken to steady-periodic state using time steps equal to 2.5×10^{−4} s (for 10, 25 and 40 Hz), 1.25×10^{−4} s (for 60, 80 and 100 Hz) and 6.25×10^{−5} s (for 150 and 200 Hz). These time steps are selected in order to accurately capture the sinusoidal drag forces in the various frequency ranges. Owing to the low Reynolds numbers of the flows calculated, the steady-periodic condition is established within a few cycles as verified by time records of the drag force acting on a spherically shaped neuromast.

A grid refinement analysis was carefully conducted prior to performing the bulk of the calculations. For this a canonical case was chosen corresponding to a tube segment with *L*_{t}/*D*=8 containing a sphere with *d*/*D*=0.50 in an oscillating flow with =4. The calculation grid was constructed through trial and error to best resolve the sheared flows near all solid surfaces. The finest grid explored consisted of 418 944 hexahedral nodes. A plot of the drag force acting on the sphere shows (not provided) that the extrapolated grid-independent result can be obtained to within 2.5 per cent on a grid with 92 728 hexahedral nodes. All other calculations were performed with this number of grid nodes.

The results of the model evaluation are shown in figure 4 for a tube-shaped canal with *L*_{t}/*D*=16 and in figure 5 for one with *L*_{t}/*D*=8. The N–S drag force solutions in figure 4 are from Part I and correspond to two neuromasts shaped as parallelepipeds with dimensions 150×150×150 μm^{3} and 100×50×50 μm^{3}, respectively, in a canal with *D*=0.25 mm. The corresponding model results are obtained by approximating the parallelepipeds as spheres with *d*/*D*=0.75 and 0.31 such that their volumes and surfaces closely correspond to those of the parallelepipeds. The N–S drag force solutions shown in figure 5 were obtained in the present study and correspond to spherically shaped neuromasts with *d*/*D*=0.25, 0.50 and 0.75. As explained in the derivation of the analytical model, the calculation of the drag force acting on a spherically shaped neuromast in a tube-shaped canal uses the concept of a hydraulic diameter to characterize the dimension of the space between the sphere and the tube. Such an approach ignores exactly where the sphere is located relative to the tube axis and two extreme possibilities arise giving the same hydraulic diameter: (i) the sphere is located in the centre of the tube, and (ii) the sphere is located on the tube wall. To check the effect this has on the numerically calculated drag force, two N–S solutions are obtained for the case *d*/*D*=0.50: one with this sphere in the centre of the tube and one with this sphere on the tube wall (figure 5). The result is a 21.4 per cent reduction in the drag force acting on the sphere located on the wall relative to the sphere located in the centre of the tube when averaged over the range of examined. This drag force reduction is due to the lower overall shear experienced by the sphere located on the tube wall (owing to its closer proximity to and interaction with low speed fluid near the tube wall) relative to that experienced by the sphere in the centre of the tube. Owing to its reliance on the concepts of a hydraulic diameter and a single characteristic velocity to determine the drag force, such an effect cannot be distinguished by the analytical model.

The results in figure 4 show fairly good agreement between model and N–S calculations of the drag force, the average differences over the range of examined being 8 and 22.9 per cent for the small (*d*/*D*=0.31) and large (*d*/*D*=0.75) neuromast, respectively. Similarly, the average differences between model and N–S calculations shown in figure 5 are 14.2, 30.8 and 65.7 per cent for *d*/*D*=0.25, 0.50 and 0.75, respectively, when the spheres are located on the tube axis. Based on the N–S case solved for a sphere with *d*/*D*=0.50 located on the tube wall, the corresponding per cent differences between model and N–S calculations for the three *d*/*D* cases with the spheres located on the tube wall are calculated to be 5.6 per cent (estimate), 12.1 per cent (exact) and 25.8 per cent (estimate). The better agreements obtained between model and N–S calculations, when the spheres are located on the tube wall for the latter, point to the goodness of the characteristic velocities (equations (2.14)–(2.16)) used to estimate the drag force in the model.

Comparisons between model and N–S drag force calculations of spherically shaped neuromasts in a tube-shaped canal with *L*_{t}/*D*=16 (not provided here) show the same trends and yield essentially the same average per cent differences as those obtained for *L*_{t}/*D*=8.

The conclusion is that for all values of the parameters *L*_{t}/*D* and explored with *d*/*D*≤0.50 the analytical model yields results in very good *quantitative* agreement with the N–S solutions. In addition, for all values of the parameters *d*/*D*, *L*_{t}/*D* and the *qualitative* agreement between model and N–S solutions is very good, with a tendency for the model to underpredict the true value of the drag force for values *d*/*D*>0.50. This is a very gratifying outcome given the approximations made in deriving the model, especially since the model allows calculations in a few seconds that would take several hours for the N–S solver.

## 5. Results

In §4.4, we established the accuracy of the physical model by reference to numerical solutions of the N–S equations for the oscillating flow in a tube-shaped canal containing neuromasts approximated as parallelepipeds and spheres (figures 4 and 5). We now use the validated model to explore the parameter dependencies of the dimensionless drag force given by equations (3.10)–(3.12) that are all of the form . For both the analytical model and the N–S numerical calculations of this study, we have varied the kinetic Reynolds number over the range 1≤≤20. Using the physical properties of water at 21°C (density, *ρ*=997 kg m^{−3}; dynamic viscosity, *μ*=9.8 10^{−4} kg m^{−1} s^{−1}; kinematic viscosity, *ν*≡*μ*/*ρ*=9.83 10^{−7} m^{2} s^{−1}) this corresponds to a frequency range 10≤*f*≤200 Hz in a canal with *D*=0.25 mm. For the analytical calculation values of *L*_{t}/*D*=4, 8, 16 and 32 are explored for 0≤*d*/*D*≤1. For the numerical calculation values of *L*_{t}/*D*=8 and 16 are explored for *d*/*D*=0.25, 0.50 and 0.75.

(<4 in sub-segment *L* and <4 in sub-segment *d*). In this case equation (3.10) applies, and the calculations of the dimensionless drag force are provided in figures 6 and 7. The plots in figure 6 show the drag force as a function of dimensionless neuromast size, *d*/*D*, in canals of different dimensionless length, *L*_{t}/*D*, ranging from 4 to 32 for a flow with =1. These plots show that decreasing the dimensionless canal length works to increase the drag force while simultaneously displacing the maximum value of this force towards smaller values of *d*/*D*. The plots in figure 7 show the variation of the drag force with *d*/*D* for different values of ranging from 1 to 4 in a canal with *L*_{t}/*D*=8. For a fixed value of *d*/*D* the drag force is seen to increase slightly with increasing kinetic Reynolds number, and the force maxima peak in a very narrow range of neuromast size centred at *d*/*D*≈0.70 for the range of explored.

(>4 in sub-segment *L* and <4 in sub-segment *d*) and case 3.2 (>4 in sub-segment *L* and >4 in sub-segment *d*). Figure 8 shows the plots of the drag force as a function of *d*/*D* for different values of *L*_{t}/*D* for conditions corresponding to case 3.3 with =8. (Results for case 3.2 are also shown but, as explained earlier, these coincide with those for case 3.3 in the range of *d*/*D* for which equation (3.11) is valid.) These plots display the same tube length dependence as in figure 6 for case 3.1. Figure 9 shows the drag force plotted versus *d*/*D* for different values of in a canal with *L*_{t}/*D*=8. These plots show that there is a critical value of *d*/*D*≈0.78 below which the drag force decreases, and above which it increases, with increasing kinetic Reynolds number. Also noteworthy is the shift of the force maxima towards larger values of *d*/*D* with increasing kinetic Reynolds number.

### 5.1 Variation of drag force amplitude and phase angle with kinetic Reynolds number

Model predictions of the dimensionless drag force plotted as a function of for different values of *d*/*D* in a tube with *L*_{t}/*D*=8 are shown in figure 10. In these plots, equation (3.10) is used for values ≤4 and equation (3.12) for values >4. For all the *d*/*D* cases examined, when ≤4 the drag force reveals a relatively weak dependence on the kinetic Reynolds number. However, when >4 and *d*/*D*≤0.75, the drag force is seen to decrease with increasing kinetic Reynolds number and the neuromast–canal system behaves as a low-pass filter. Over the entire range of explored, the drag force increases with increasing *d*/*D* up to *d*/*D*=0.8, approximately, as of which the drag force decreases with further increases of *d*/*D*. This behaviour is in accordance with the variations displayed in figures 7 and 9.

Model predictions of the drag force phase angle obtained using equation (3.14*b*) are shown in figure 11. For fixed *d*/*D* the phase angle increases with increasing kinetic Reynolds number, from zero at =0 to a value that would be 90° for →∞. For fixed kinetic Reynolds number the phase angle increases to an asymptotic value that is attained by *d*/*D*=0.8, approximately.

## 6. Discussion

In §1 of this study we asked two interrelated questions. Let us now consider them in light of the insight provided by the analytical model.

How does the drag force acting on a neuromast inside the LLTC vary as a function of the parameters affecting it? In particular, is there a purely physical optimal neuromast size? The physical–mathematical analysis shows that three dimensionless groups (

*d*/*D*,*L*_{t}/*D*, ) independently affect the dimensionless drag force acting on a neuromast inside the LLTC of a fish (see equation (3.14*a*)) and the model calculations reveal the following noteworthy trends.For viscous-dominated flows inside all canal sub-segments (case 3.1), increasing the total length of the canal decreases the drag force and simultaneously displaces the force maximum towards larger values of

*d*/*D*(figure 6). Both of these effects are due to the increased pressure drops in the canal sub-segments of length*L*relative to that in the sub-segment of length*d*. Also, for a fixed total canal length,*L*_{t}, increasing the kinetic Reynolds number, , only slightly increases the drag force without significantly affecting the*d*/*D*location of the force maxima (figure 7).For inertia-dominated flows inside the canal sub-segments of length

*L*and inertia or viscous-dominated flows inside the canal sub-segment of length*d*(cases 3.2 and 3.3), the same dependence on the total canal length,*L*_{t}, as discussed above for case 3.1 applies (figure 8). However, when the total canal length is fixed and the kinetic Reynolds number is varied (figure 9), two noteworthy differences arise relative to case 3.1: (1) increasing the kinetic Reynolds number significantly displaces the force maxima towards larger values of*d*/*D*, and (2) there is a critical value of*d*/*D*≈0.78 below which the drag force decreases, and above which it increases, with increasing kinetic Reynolds number.In considering the question of the existence of a purely physical optimal neuromast size, we exclude (undeniably important) biological considerations in order to focus on the purely physical part of the question since, irrespective of the biology, the sensory organ must obey the physics. The data presented in figures 6–9 illustrate an important point for all the cases examined. This is that the drag force acting on a neuromast in a canal must decrease towards zero for small and large values of

*d*/*D*. For small*d*/*D*this is because the neuromast surface and volume are proportional to*d*^{2}and*d*^{3}, respectively, and as*d*decreases so must the friction and inertia drag forces acting on the neuromast; see equation (3.3). For large*d*/*D*this is because the velocity and acceleration of the flow through the canal, on which the same friction and inertia forces depend, are significantly reduced by the large pressure drop caused by the neuromast blocking the flow; again, see equation (3.3). As a consequence, the drag force acting on a neuromast must, indeed, display a purely physical maximum which for the conditions of this study occurs in the range 0.65≤*d*/*D*≤0.85, depending on the nature of the flow in the tube sub-segments. The*d*/*D*locations of the force maxima vary with both*L*_{t}/*D*and for inertia-dominated flows (cases 3.2 and 3.3) but primarily with*L*_{t}/*D*for viscous-dominated flows (case 3.1).The kinetic Reynolds number, , is a dimensionless flow oscillation frequency. Figure 10 shows that, irrespective of neuromast size

*d*/*D*, for values <4 the neuromast–canal system is relatively insensitive to the value of frequency. By contrast, for values >4 and neuromast sizes*d*/*D*<0.75, approximately, the system behaves as a low pass frequency filter, the effect being most noticeable for small values of*d*/*D*.By its definition, equation (3.9

*b*), the drag force phase angle,*β*, is a measure of the effects of inertia forces relative to viscous forces acting on a neuromast, small values of*β*corresponding to viscous-dominated drag forces and large ones to inertia-dominated drag forces. Figure 11 shows that relatively small increases in either the neuromast size,*d*/*D*, or the flow kinetic Reynolds number, , lead to significant increases in the inertia drag forces associated with flow accelerations. For a phase angle corresponding to*β*=45° the quantity*Bω*/*A*in equation (3.9*b*) is equal to unity, denoting equal contributions of inertia and viscous effects to the drag force.

What are some of the possible implications to biology, and to the design of corresponding bioinspired neuromast-like flow motion sensors, of the drag force parameter dependencies observed?

Whether or not nature takes advantage of the fact that the dimensionless drag force acting on a neuromast maximizes at particular values of

*d*/*D*(depending on*L*_{t}/*D*and ) for the moment remains unknown. (To resolve this issue would require the statistical examination of a large, representative body of data pertaining to fish neuromast and canal dimensions and the associated flow oscillation frequencies and pressure gradients they are exposed to.) The fact that this ‘physically optimal’ size exists theoretically does not necessarily mean that neuromast–canal systems have evolved towards it in practice because other important factors are simultaneously at play. For example, they may relate to the material properties, structure, and shape of the neuromast, which may endow it with low threshold sensitivities relative to the applied drag force. Consequently, the neuromast may not need to be physically optimal in terms of its size to maximize drag. The story is different, however, for a synthetic sensory organ modelled after the LLTC. In this case, the sensor (‘neuromast’) material properties, structure and shape are bound to be vastly simpler relative to the biological, and it is both advantageous and possible to capitalize on designs that emphasize a range of*d*/*D*maximizing the drag force in the range of kinetic Reynolds numbers of interest.In relation to the above point, van Netten (2006) has elegantly shown that the determination of the optimal frequency selectivity of a neuromast requires an analysis of the mechanics of neuromast cupula excitation based on the equation of motion for the displacement of the cupula. Such an equation includes a drag force term acting on the cupula due to the relative motion between the fluid and the cupula. van Netten (2006) reports that the cupula displacements used in his analysis were induced by vibrating a sphere placed a few millimetres from the cupula in the supra-orbital canals of the ruffe (

*Acerina cernua*or*Gymnocephalus cernuus*) and the African knife fish (*Xenomystus nigri*). Such a study is very valuable but, because it does not resolve the drag force acting on the cupula*as a function of the canal parameters d*/*D*,*L*_{t}/*D*and , it provides limited insight concerning the dependence of neuromast sensory performance on*these*interrelated parameters. To elucidate these dependencies requires repeating van Netten's analysis using equation (3.8) of the present study to evaluate the time-varying drag force. However, such an effort is beyond the scope of the present investigation and remains to be done.Since in this study, the value of

*d*maximizing the drag force acting on a neuromast can be obtained though different (*h*_{N},*d*_{N}) combinations corresponding to a half-prolate spheroid, or through different (*h*_{N},*w*_{N},*l*_{N}) combinations corresponding to a more general half-spheroid. Clearly, these degrees of freedom are advantageous to the evolutionary adaptation of the sensory organ which must optimize many other factors, related and not related to the drag force acting on it, that affect its overall performance characteristics. Similarly, since , different (*ω*,*D*) pairs can give the same value of . For example, for the case of water oscillating with a value of =2 in a canal with*D*=0.5 mm the flow oscillation frequency is*f*=5 Hz, while in a canal with*D*=0.25 mm the flow oscillation frequency is 20 Hz. Notwithstanding, because the results of this study are provided in dimensionless form, they apply to a large range of geometrical and dynamical conditions. For the same reason, any variations of kinematic viscosity, whether intrinsic to the canal fluid or due to variations with temperature (Esther*et al*. 2000), are captured by the dimensionless kinetic Reynolds number, .The model equations for the drag force acting on a spherically shaped neuromast are of the form given by equation (3.14

*a*). Whether applied to a biological or a synthetic system, they provide a quick and relatively accurate way to calculate relative parameter dependency effects. In this respect, the equations are especially valuable for synthetic sensor design. For example, for the dimensionless canal lengths,*L*_{t}/*D*, and kinetic Reynolds numbers, , examined in this study, the dimensionless drag force maximizes for values of*d*/*D*ranging between 0.65 and 0.85. Clearly, for these ranges of*L*_{t}/*D*and , it would be advantageous to design a synthetic neuromast with 0.65<*d*/*D*<0.85, approximately, for maximum physical performance. In particular, for flows with kinetic Reynolds numbers >4, a neuromast with*d*/*D*≈0.78 will experience a drag force that is essentially*independent*of the kinetic Reynolds number (figure 9) and of magnitude very close to its maximum value.

## Acknowledgments

Much of the material contained in this paper was developed while the author was a Guest Professor in the Department of Neurobiology and Cognition Research at the University of Vienna during the summers of 2007 and 2008. The author expresses his sincere appreciation to Professor F. Barth for this courtesy appointment and for many exciting and illuminating discussions on animal sensory systems with him and his students. Many thanks go to Dr Luis Rosales for his assistance with the N–S based numerical solutions of the flow oscillating about a sphere inside a tube. The study was performed under a DARPA BioSenSE Award via AFOSR grant no. FA9550-05-1-0459.

## Footnotes

↵The drag force of interest acts on the neuromast cupula, that part of the neuromast exposed to the pressure and shear forces induced by the flow past it. In this work, we shall simply refer to the ‘neuromast’, having removed any ambiguity concerning where the drag force is applied.

↵This check is readily performed by using equations (2.14)–(2.16) to estimate

*V*_{S}, which is equivalent to assuming the sphere is fixed in a fluid with oscillating velocity*V*_{d}.↵Because these specifications for

*δ*are restricted to the*d*/*δ*terms in equation (3.3), the corresponding specifications for in equation (3.13) apply only to the underlined terms of that equation.↵If instead of the cosinusoidal form of equation (3.4) a sinusoidal form is used, equation (3.8) becomes

*F*_{S}=|*F*_{S}|sin(ω*t*+β) but equations (3.9*a*) and (3.9*b*) remain unchanged.↵See footnote 3 for a clarification of the underlined terms in equation (3.13).

↵In this subsection, we depart from the kinetic Reynolds number notation noting that |

*V*_{L}|≈*ωD*and, therefore, that .- Received July 9, 2008.
- Accepted September 19, 2008.

- © 2008 The Royal Society