Bio-inspired all-optical artificial neuromast for 2D flow sensing

We present the design, fabrication and testing of a novel all-optical 2D flow velocity sensor, inspired by a fish lateral line neuromast. This artificial neuromast consists of optical fibres inscribed with Bragg gratings supporting a fluid force recipient sphere. Its dynamic response is modelled based on the Stokes solution for unsteady flow around a sphere and found to agree with experimental results. Tuneable mechanical resonance is predicted, allowing a deconvolution scheme to accurately retrieve fluid flow speed and direction from sensor readings. The optical artificial neuromast achieves a low frequency threshold flow sensing of 5 mm s−1 and 5 μm s−1 at resonance, with a typical linear dynamic range of 38 dB at 100 Hz sampling. Furthermore, the optical artificial neuromast is shown to determine flow direction within a few degrees.


Introduction
The design of our optical sensor is directly inspired by a sensory modality used for flow detection by aquatic vertebrates. Fish and amphibians have an array of discrete mechanical sensors at their disposal called neuromasts, which are distributed along the head and trunk. With these neuromasts they can perceive a 1D projection of local fluid motion or flow relative to the body (Dijkgraaf 1963). The biophysical properties of neuromasts, such as the interplay of physiology, mechanics and fluid dynamics, have been studied intensively elsewhere (van Netten 2006, van Netten andMcHenry 2013).
There are two types of neuromasts, each with their beneficial physical properties to help the fish perceive freestream (DC) and dynamic or oscillatory (AC) flow (Bleckmann et al 2004). Superficial neuromasts (SN) are present on the surface of the body and are in direct contact with the surrounding medium. They are tailored to perceive steady (DC) and low frequency fluid flow velocity. The canal neuromasts (CN) are not in direct contact with the freestream flow, but are housed in internal canals. They are deflected through the pres sure difference via flexible membranes or pores in these canals, thereby effectively perceiving freestream (AC) fluid acceleration.
The perceived local fluid flow at each neuromast over time can be concatenated to a spatiotemporal flow pattern, which augments the fish sensory perception (Curcic-Blake and van Netten 2006). This enables fish to sense fluid perturbations generated by moving sources. In fish behavioural experiments, the lateral line has been shown to be instrumental in many specific behaviours, for instance, prey detection, schooling behaviour, and spatial orientation (Coombs and Montgomery 1999, Ghysen and Dambly-Chaudiere 2007, Tao and Yu 2012. In order to mimic this biological near-field sensing, several implementations of artificial neuromast sensors have been developed, which also measure a 1D projection of local fluid flow. Some sensors make use of hot wire anemometry . Here, a suspended hot nanowire is cooled down by fluid flow. This links a measureable change in the temperature dependent resistance to fluid flow speed. Most artificial SNs, however, rely on sensing generated strain at the base of a deflecting lamella or canti lever structure in response to fluid flow. This design is favoured since the protruding structures escape unwanted boundary layer effects. Several techniques used for strain sensing include lamella mechanical micro-sensors (MEMs) (Fan et al 2002, McConney et al 2009, Yang et al 2007, 2010, ionic polymer-metal composites (IPMC) (Abdulsadda and Tan 2012, Chen et al 2013 and soft polymer membranes without (Asadnia et al 2013) and with cantilever structures (Kottapalli et al 2014. A technical review of recent contributions to the field of artificial neuromast and artificial lateral lines (ALL) shows an increasing interest in this field (Liu et al 2016). In some cases, orientations of 1D-sensitive sensors are alternated to sense multiple projections, allowing measuring flow perpendicular to the array, e.g. Yang et al (2010) and Ahrari et al (2017).
These SN designs rely on electric methods to operate in wet conditions. Because electric signals are susceptible to noise pickup over large distances, deployment of remote and large scale artificial lateral lines (ALL) has been somewhat limited. We therefore present an all-optical artificial neuromast which aims to address these issues. Our design utilises fibre Bragg gratings (FBGs), which allows the sensor data to be transmitted through fibre optic cables, thereby enhancing the scalability for all-optical ALLs. Furthermore, our design enables two dimensional fluid flow measurements, thus increasing the information per point measurement.
FBGs are sections in an optical fibre that have been modified to include a periodic variation in the refractive index along the fibre length. This structure reflects light at a specific wavelength that is determined by the spacing of the FBG structure. Stretching or compressing the fibre changes the FBG structure spacing and therefore increases or decreases the reflected wavelength. If the FBG is illuminated with a broadband optical source containing a wide range of wavelengths then the reflected (Bragg) wavelength can be interpreted as a function of the applied strain.
Different geometries consisting of two or more optical fibres glued together have been used in order to determine the curvature of the end position of the combined fibre structure (Araujo et al 2002). This requires the gratings to be located outside the neutral (bending) axis of a cantilever structure. Examples of FBG curvature sensors (Flockhart et al 2003) and accelerometers (Fender et al 2007) have also been developed using multicore optical fibres with multiple cores positioned away from the neutral axis.
Cantilever Bragg grating flow sensing has been used for monitoring steady flow rates (Lu and Chen 2008) and flow perturbations in response to a bluff body (Takashima et al 2004) in pipes. However, the dynamic properties of the sensors were not examined.
We first present a combined hydrodynamics and strain-structure model which enables sensor characteristics such as its mechanical sensitivity and frequency response to be predicted. Using steady state contact deflection, we infer the linear dynamic range of the sensor. Finally, through hydrodynamically stimulating the sensor at different frequencies, we verify the sensor characteristics and employ a related method to reconstruct flow speeds from sensor readings.

Sensor design
The sensor physically resembles a fluid force recipient spherical body and a fibre support structure providing elastic coupling (figure 1(a)). To model the signal of a deflecting sensor, the elastic support is treated as an endloaded cantilever beam. Using Bernoulli's beam equations, we model the support with length h, having a circular cross section, and possessing flexural stiffness EI (with Young's modulus E and the second moment of area I).

FBG sensing
Equations (2.1) and (2.2) show the relationship between the force F applied at a beams tip, the resulting tip displacement magnitude ∆τ and the generated strain ε at a distance d from the neutral bending axis at a height z from the fixed cantilever end (Benham et al 1996). The beam is compressed in the bending direction, which decreases strain. On the opposite side of the bending axis, where the sign of d is opposite, the beam stretches locally, increasing strain. The all-optical sensor has a height of 64.8 mm and spherical body radius of 4.00 mm. The fibre structure is glued into a 3D printed block, allowing four separate outgoing optical fibres.

∆τ =
This strain is measured locally via FBGs (figure 1(a)). The FBG comprises a periodic refractive index modulation with a period on the scale of the wavelength of light, such that it reflects light at the Bragg wavelength, λ B , which is twice the inter-grating distance or grating period Λ.
The Bragg wavelength changes as a function of strain (equation (2.3)) and temperature (Hill and Meltz 1997). The thermal effect can be compensated using differential strain configurations in which the thermal effect is common to all considered gratings (Flockhart et al 2003). Both the local Bragg wavelength shifts (i.e. sensor signal) and sensor deflections can be practically obtained to ascertain to what extent these are consistent with our present model. Common FBG wavelengths are in the communications band and we choose to work with gratings in this region of λ B = 1560 ± 40 nm.

Isotropic sensor mechanics
The support comprises four discrete standard communication SMF-28 fibres, where the light guiding core is centred in a silica outer cladding. The mechanical properties of the fibre support structure are largely determined by the bending stiffness K which is, apart from sensor height, governed by the geometry and composition of the cross section of the elastic support.
Since the support matrix and fibres form a composite material, its flexural stiffness EI depends on the mechanical properties of both materials. The added flexural stiffness of each fibre depends on its squared distance to the neutral bending axis (Benham et al 1996). Although the bending axis may rotate by deflecting the sensor in different directions, we show below that the flexural stiffness-and therefore the bending stiffness K-is independent of a variable bending axis offset α.
In order to show that the flexural stiffness is circular symmetric, only fibres 1 and 2 are considered, since the cross section is mechanically symmetrical above and below the bending axis. When the sensor is deflected over an axis with a variable offset angle α compared to the square lattice, the individual distances d 1 , d 2 between the fibres and the bending axis change. But since they form sides of an identical right angle triangle, their summed squared distance to the bending axis (d 2 1 + d 2 2 = p 2 /2) remains constant and is independent of α. It is only dependent on a fixed p, the distance between two adjacent fibre cores. The flexural stiffness of the support is therefore mechanically circularly symmetric; no preferred bending direction exists.
This results in a direction independent, or isotropic, bending stiffness K determined at the tip of the fibre structure where r c and r s denote the cladding and support radii and E c and E s their respective Young's moduli.
Similarly, in this square lattice configuration (figure 2(b)), we show below that the magnitude of a differential strain vector − → δ , generated by a fixed magnitude of tip deflection ∆τ , is not affected by a variable bending axis offset α. First, we can define the fibre distances d 1 = cos (π/4 − α) When deflecting the sensor, with an offset bending axis α, in the direction α + π/2 and using equation (2.2), the Cartesian differential strain projections then become . Four fibres are embedded in an adhesive matrix. In the two fibres on the right, the core is made visible by propagating white light in the fibres to illuminate the core.
where n denotes the generated strain at core number n. Then, working out the magnitude and angle of differential strain in polar coordinates, we find that (2.9) The absolute differential strain magnitude r δ shows that we can generalise the use of individual fibre distances d to the inter-fibre distance p in case of differential strain c.f. equation (2.2) & (2.8). Furthermore, the direction θ δ of increased differential strain lags the bending axis by π/2 and is therefore opposite to the bending direction, which is consistent with equation (2.2).
Since the strain vector − → δ is opposite to to the bending vector (figure 2(b)), and we aim to measure flow speeds in the bending direction, we choose a pairing of sensor cores such that the projection of differential wavelength shift ( − → δλ ∝ − − → δ ) is in the bending direction. With respect to the sensor orientation as depicted in figure 2, the Cartesian projections of wavelength difference are given by where c n denotes the generated Bragg wavelength shift at core number n. Alternative combinations of cores, such as along the diagonals of the lattice, also provide the required orthogonality. This diagonal combination requires four functional cores and is therefore less redundant.
The induced differential wavelength shift − → δλ and deflection − → ∆τ can conveniently be defined as vectors and are thus only a dimensionless linear geometrical factor apart (equation (2.12)). (2.12) Using practical values for p of tenths of millimetres, combined with sensor height h in the order of centimetres, the value of G is typically in the order of 10 −7 .

Dynamic properties
The physical parameters of both the elastic support and sensor body can be adjusted to predict and obtain desirable dynamic and filter characteristics. We obtain the frequency dependent sensitivity of the sensor by adapting a model for CNs found in fish (van Netten 2006). The magnitude of the complex frequency response FR( f ) is defined as the frequency dependent ratio of sensor response (in this case sensor body motion) per unit fluid velocity. Its argument constitutes the frequency dependent phase lag. In this model, given the relative dimensions of the sensor body and fibre support, we neglect fluid forces acting on the support.
Only three independent physical parameters determine the frequency response: the transition frequency f t , which selectively shifts a constant shaped Bode plot along the frequency axis, a resonance number N r , which determines the shape of the associated Bode plot, and the novel buoyancy factor b of the sensor body, which also affects the Bode plot shape. The frequency response and its parameters are described by equations (2.13)-(2.16). (2.14) (2.16) The sensor is a dampened resonator; its resonant character is parameterized by N r and affected by the bending stiffness K, sensor body radius a, fluid viscosity μ, fluid density ρ fluid , and body density ρ body .
The resonance frequency f r (equation (2.17)) indicates the frequency at which the maximum of the FR is located. The resonating behaviour is quantified through the quality (Q) factor (equation (2.18)), which only depends on N r . Keeping the Q factor low causes the sensor low pass filter function (i.e. frequency response) to remain as flat as possible, effectively increasing its useable bandwidth. (2.18) The bandwidth itself is bounded by a cut-off frequency f c , the frequency at which the response in the fall-off region matches the DC response and is given by . (2.19)

Design optimisation
Using the dynamic and mechanical properties, we can optimise some sensor parameters with respect to fluid sensing. With FR( f = 0), the response at low frequencies can be found. Together with the cut-off frequency this leads to a sensitivity-bandwidth (SB) product, which only depends on b: . (2.20) The SB is clearly maximal when b = 0, then SB = √ 3/(2π), which is an 73% increase compared to neutrally buoyant neuromasts as found in fish (van Netten and McHenry 2013). Using polyethylene spheres with an effective b of about 0.05 in water at room temperature, allows for an increase of SB of about 65% over neutrally buoyant sensors.
Given the geometric factor (equation (2.12)) and the frequency response (equation (2.17)) at low frequencies, we can also optimise the generated differential wavelength shift per fluid velocity by varying the inter-fibre distance p.
In equation (2.21), it is reflected in the numerator that the wavelength shift per fluid velocity increases linearly with the core distance p. However, the optical fibres have to be compressed and stretched at effectively larger distances, adding to the flexural and therefore bending stiffness which is reflected in the denominator with a factor of p 2 . As a consequence, p has an optimal value for maximising detected wavelength shift per fluid velocity at With the chosen materials and related Young's moduli, this optimal distance is smaller than the core diameter. Therefore the optimal design has the fibres as close as possible to each other. Fabrication constraints limited the practical separation of the fibres to a minimum of p = 0.3 mm.

Methods
In the manufacturing process, four standard communication optical SMF-28 fibres (jackets removed, r c = 62.5 μm, E c = 75 GPa) are suspended in a square column formation. This leaves room for the optical adhesive to flow around the fibres. A custom glue dispenser and UV-curing system is encapsulating the four fibres and is slowly moved along the fibre formation. This embeds the fibres in a matrix (E s = 0.14 GPa) of UV-cured optical adhesive (NOA68, Norland Products Cranbury, NJ, USA). A cross section of the resulting fibre structure is shown in figure 2(c). The fibre structure is then glued into a mounting block using an epoxy adhesive, in which a small diameter hole acts as a fixed cantilever point for the sensor. The sensor with length h = 64.8 mm is fitted with a polyethylene sphere a = 4.00 mm, b = 0.05. Figure 2(c) shows a cross section of the tip of the support structure, which is the closest indication of the cross section at the FBG height z. Here, the average distance between the fibre cores is found to be p = 313.1 ± 5.4 μm. The optical fibres do not form a perfect square, therefore we can expect some variation of bending stiffness depending on the bending direction. With the average radius of the support r s = 277.4 ± 5.5 μm, and the centre of the inscribed FBG sections located at z = 5 mm, we expect a geometrical factor (equation (2.12)) of G = 3.22 ± 0.05 · 10 −7 , which describes the relation between the generated differential wavelength shift (sensor signal) per deflection at height h. On the basis of the dynamic properties described in section 2.4, we expect the sensor to have a resonance frequency of 13.4 Hz and a Q factor of 13.4.
First, we validated that the sensor response is linear in our operational range by mechanically deflecting the sensor. This allows for a linear relation between sensor body motion and resulting optical signals via a geometric factor, and is a requirement for acquiring the frequency response FR using the current method. Using a linear stage, the sensor was deflected in steps of 0.5 mm while monitoring the sensor signal. Then, in order to infer its dynamic properties, i.e. the FR, the sensor motion and sensor signal were monitored in response to a hydrodynamic dipole stimulus (figure 3).
Sensor body motion was measured with a calibrated Zeiss Axiotron microscope. A 2D position sensitive detector (On-Trak PSM 2-2) allowed for high-speed 2D tracking of an attached reflectance marker within a range of 500 μm × 500 μm in the objective focal plane. Strain-induced FBG wavelength peak shift was measured with picometer resolution at 5 kHz using an optical interrogator (Micron Optics si225, Atlanta, USA). The difference in reflectance peaks of cores 1 and 2 and those of core 4 and 1 (figure 2(a)) are used as Cartesian x and y projections of measured wavelength shifts (equations (2.10) and (2.11)).
A hydrodynamic dipole source was produced using a Bruel & Kjaer 4810 mini-shaker driving a submerged sphere ( = 9.9 mm). This stimulus was levelled with the sensor height h. The water tank setup (figure 3) was placed on a vibration isolation table (Newport VW series) as to avoid mechanical noise. Sinusoidal stimuli were generated, and sensor body motion synchronously acquired, using a CED power1401 data acquisition system. The mini-shaker was calibrated for frequencies ranging from 1 to 200 Hz, with a constant travel amplitude of about 58 μm. A model for viscous flow (van Netten 2006) was used for calculating the fluid flow velocity produced by the calibrated stimulus. The stimulus sphere was positioned at a distance r, measured using a micrometer stage.
The sensor was pre-stimulated diagonally at a given stimulation frequency for at least 3 seconds before sampling, in order to avoid transients. For both sensor motion and sensor signal, we obtained the amplitude and phase lag by applying a flat-top window suitable for low amplitude and resolution data (Heinzel et al 2002), and calculating the discrete Fourier transform. From the resulting spectrum magnitude, we take the maximal value at the stimulation frequency as the response amplitude. The corresponding imaginary part yields the phase-lag at that frequency. In order to prevent spectral aliasing in Fourier analysis, 64 low pass filtered periods of sensor motion with 512 samples per period were sampled for a given stimulation frequency.

Results
The Cartesian (x, y) sensor signals for a sensor at rest showed normal distributed noise on the measured differential wavelengths with a standard deviation of 2.02 pm at 5 kHz sampling, which defines a lower bound for the dynamic range of the sensor.

Linearity
To test for linearity, the sensor has been deflected in steps of 0.5 mm up to 10 mm (black) and back to origin (cyan) for both positive and negative deflections, where each position was held for at least four seconds.
From figure 4, we infer that the sensor response is linear for deflections up to 5 mm as a conservative  Bioinspir. Biomim. 13 (2018) 026013 estimate (R 2 > 0.999), and approximately linear up to 10 mm with some small deviations. On average, a deflection of 5 mm corresponds to a differential wavelength shift of 1.78 nm. The slope of the response matching this linear part yields a measured geometric factor (equation (2.12)) of 3.56 · 10 −7 under static deflection conditions.

Dynamic response
To determine the frequency response and further verify the geometric factor at dynamic submillimetre deflections, the sensor was hydrodynamically stimulated at frequencies from 2 to 100 Hz. Figure 5 shows the frequency response for the sensor body motion in the x and y directions.
Here, the x (black) and y (cyan) projections of measured sensor body motion (squares) and phase (triangles) are plotted for each stimulation frequency. Their respective fits result in the fitted param eters b = 0.05, f t = 0.021, N r,x = 7.12 · 10 4 and N r,y = 6.52 · 10 4 .
The measured frequency response (figure 5) shows a slight difference in resonance peaks and phase lag for the x and y directions. From this measurement, we find that the sensor has a resonance frequency of 14.7 Hz in its x direction, with a Q factor of 15.3. In the y direction, we find a resonance frequency of 14.0 Hz and a slightly lower Q factor of 15.1. This slight difference is most visible near the phase flip and amplitude maxima and is reflected in their different fitted values for N r . This difference can therefore be attributed to the mechanical properties of the non-perfect square lattice inside the fibre support structure; a slight directional difference exists. Figure 6 shows the wavelength shift signal and body motion over time of five trials of the frequency response measurement. As expected, the relative orientation and magnitude of both types of data match up. From the ratio of sensor signal amplitude to sensor motion amplitude, we obtain a measured geometric factor of 3.58 · 10 −7 in the submillimetre domain under dynamic conditions.
Due to the small directional difference in bending stiffness and resulting frequency response, the x and y projections of motion, and therefore sensor signal, have a slight phase lag difference. This produces sensor motion in ellipses rather than straight lines at stimulation close to either x or y resonance frequencies. In the next section we demonstrate a method to correct for this phenomenon.

Discussion
The measured sensor behaviour has been shown to compare well with the modelled input-output and frequency characteristics. The sensor acts as a dampened resonator both driven and dampened by fluid forces. The results show that the measured frequency response (figure 5) is accurately described by the combined hydrodynamics and strainstructure model. Both large, static, mechanical deflection (G = 3.56 · 10 −7 ) and small, dynamic, hydrodynamical stimulation (G = 3.58 · 10 −7 ) aligns with the modelled mechanical properties of the sensor based on the cross section from the sensor tip (G = 3.22 ± 0.05 · 10 −7 ).
This model can therefore be used to relate sensor signal via sensor body motion to local fluid flow in two steps. First, the sensor signal can be translated to sensor body motion with the geometric factor G. Then, by a deconvolution via the frequency response (equation (2.13)), which incorporates the fluid-dynamic characteristics of the sensor, we can accurately obtain the local 2D flow velocity (see figure 7).

Sensitivity and dynamic range
The measurement of the frequency characteristics in figure 5 shows an expected peak, which indicates the resonance frequency, where the sensor is most sensitive. Higher frequencies will be filtered out, i.e. there's a constant downward slope in the frequency response curve. At frequencies below resonance, the frequency response plateaus and will still pick up low frequency and DC flow (see equation (2.19) and table 1), although it is not specifically designed to do so.
The lowest detectable flow speeds of the artificial neuromasts affect the detection limits for tracking objects using artificial lateral lines (Boulogne et al 2017). We can infer this threshold velocity, at which the expected sensor signal equals noise levels, by taking the ratio of the measured noise in − → δλ to G times the frequency response (table 1).
When down sampling the signal by averaging with a factor 50, the noise levels dropped from 2.02 pm to 0.28 pm with a factor close to √ 50, as would be expected from normally distributed noise. The sensor is shown to be linear up to deflections of about 5 mm (figure 4). Taking the corresponding 1.78 nm wavelength shift and σ noise as the upper and lower boundaries respectively, this amounts to a linear Dynamic Range of 29 dB at 5 kHz sampling and 38 dB when down sampled to 100 Hz.

Two-dimensional fluid flow sensing
Two orthogonal projections of flow can be measured by the all-optical artificial neuromast using the information of at least three cores. This requires that four fibres are positioned in a square lattice embedded in a support structure. The sensor cross section as shown in figure 2(c) shows that the centres of the four fibres do not form a perfect square, so some directional variance in stiffness and therefore sensitivity exists. By employing a 2D method of sensing sensor body motion, we can measure, and correct for, this directional variance. In practice, by processing the sensor signals via a deconvolution along two orthogonal dimensions, we can translate sensor signal and body motion with their respective phase lag and resonance properties to a 2D representation of local fluid flow. Figure 7 shows the deconvolution process (for details, see Smith (1997)). Here, measured sensor signals are first band pass filtered in order to reduce high  and low frequency noise. In the frequency domain, the real sensor amplitudes are divided by the amplitudes of the FR. The FR phase is subtracted from the imaginary valued sensor phase information. Using an inverse Fourier transform, the resulting amplitude and phase information is then transformed back into the time domain. For longer or continuous sensor readings, taking a small buffer period into account, this procedure can be performed real-time using sliding time windows. In addition, we can express the uncertainty of the sensor readings both in Cartesian projections of measured fluid flow as well as a magnitude and direction representation. As is clear from figure 7(d), the reconstructed velocity contains some magnitude and direction variance. Ignoring velocity readings under 0.05 mm s −1 and taking the ratio of V x and V y over time, the variation in the reconstructed direction of velocity oscillation amounts to a standard deviation of σ = 1.8°, with an oscillating flow amplitude of 0.54 mm s −1 . The sensor is therefore able to reliably determine fluid flow direction within a few degrees.

Conclusion
We have shown that the fibre structure mechanically acts in accordance with our combined hydrodynamics and strain-structure model. The predictable nature of dynamic behaviour allows for deconvolution of measured strain induced FBG signals via the frequency dependent velocity sensitivity which, in turn, allows for direct translation to fluid velocity. The dynamic characteristics are tuneable by varying the sensor dimensions and mechanical properties, allowing a tailored design for specific use cases.
The presented strain-structure model allows optimizing the sensor design for generating strain and thus sensor signal per sensor motion. The observed geometric factor relating sensor body motion to observed differ ential wavelength shift is in accordance with the modelled geometric factor. It is consistent throughout the measurements and allows individual sensors to be calibrated by using conventional displacement sensitive methods, such as mechanical deflection or visual monitoring.
A welcome enhancement from its 1D-sensitive biological counterpart and state-of-the-art artificial neuromasts is that by interrogating the sensor in two perpendicular axes, a single sensor will provide information about both the magnitude and the angle of the local flow velocity. Although it is shown that stationary artificial lateral lines can reconstruct dipole sources and moving artificial lateral lines can detect obstacles without a second dimension, this extra information might be able to increase the ALL detection precision and effectiveness.