Floquet stability analysis of a two-layer oscillatory flow near a flexible wall

Abstract We investigate the linear Floquet stability of two fluid layers undergoing oscillations in the direction parallel to the flexible wall that separates them. This canonical configuration is inspired by the cerebrospinal fluid flow in the spinal canal of subjects with hydromyelia/syringomyelia. The analysis focuses on the marginal conditions for the onset of instability, and how these depend on the spatial wavelength of the perturbation, and on the values of the control parameters, which are the two channel widths, the Reynolds number and the wall stiffness. Unstable perturbations are found to oscillate synchronous with the base flow. The wavelength of the most unstable perturbation, of the order of the stroke length of the basic oscillatory motion, depends strongly on the wall stiffness, but is only weakly influenced by the channel widths and the Reynolds number. In general, around criticality, it was found that increasing the Reynolds number has a destabilizing effect, and that decreasing the canal widths stabilizes the instability. The wall stiffness on the other hand has a non-monotonic effect, exhibiting an intermediate value for which the instability is maximally amplified. The present analysis is a first step towards a better understanding of the physical mechanisms that govern many (bio)fluid mechanical problems that involve oscillatory flows near compliant walls.


Introduction
Understanding the role played by pulsatile or oscillatory flows in the vicinity of deformable walls is crucial in many engineering Jaffrin and Shapiro [1971] and biomedical applications Grotberg and Jensen [2004].In fact, many of the physiological systems of the human body present this kind of motion.The flow of air in the respiratory system, that of blood in the circulatory system, or the flow of cerebrospinal fluid (CSF) in the central nervous system (CNS), constitute relevant examples where the oscillating driving mechanisms, together with the interaction of the fluid with the compliant walls, determine the nature of the fluid motion Linninger et al. [2016], Heil and Hazel [2011].Nevertheless, despite their relevance, many aspects of oscillatory flow in the presence of compliant walls are still poorly understood.In the present paper we focus on the flow of CSF in the cervical spinal subarachnoid space and its interaction with the fluid in the central canal of the spinal cord through the deformation of the spinal cord tissue that separates both fluid layers (Figure 1a).We focus on this configuration in the context of two pathologies, named syringomyelia and hydromyelia.Both conditions are similar and involve the accumulation of fluid in the spinal cord.In syringomyelia it occurs in the form of a slender fluid-filled cavity, called syrinx, approximately 2-6 cm long and 3-6 mm wide (Figure 1b) Elliott et al. [2013], whereas in hydromyelia there is an abnormal dilation of the central canal that typically spans a longer section of the spine, about 3-4 vertebrae (Figure 1c) Roser et al. [2010].In hydromyelia, the dilated central canal is considered to be connected to the fourth ventricle (which is the case in newborns and young children), whereas in syringomyelia, the cavity may be isolated as a result of age-related obstruction of the central canal.Nevertheless, the distinction between both diseases is often hard to make, and some physicians might even use both terms interchangeably.The pathogenesis of these diseases is still unknown, but there is a consensus that CSF hydrodynamics may play an important role.In fact, many cases of hydro-/syringomyelia are associated with a malformation of the hindbrain named 'Type I Chiari Malformation' (CMI), in which the cerebellar tonsils descend into the entrance of the spinal subarachnoid space and partially block the natural oscillatory motion of CSF (Figure 1b).This back-and-forth motion of CSF in the spinal subarachnoid space surrounding the spinal cord is in principle coupled to the motion of fluid in the central canal, with associated periodic deformations of the compliant spinal cord tissue.Whereas in healthy subjects these deformations may remain small, in subjects with Chiari Malformation the changed hydrodynamics in the subarachnoid space may induce larger self-sustained deformations of the spinal cord tissue separating the two fluid layers.Prolonged exposure to such deformations may then in turn be associated with additional oedema or fluid transport into the fluid cavity, leading to the growth of cavities, whose size may be marked by the wavelength of the preferential deformations.As a first step to explore if this proposed mechanism is viable, in the present paper we study the underlying fluid-structure interaction problem using a simplified model description.
The main features of the anatomy of the cervical spine that are relevant to the problem at hand are given in figure 1.The spinal subarachnoid space is a thin annular canal, approximately  * SSAS ∼ 1 − 4 mm wide, bounded internally by the pia membrane (thickness  * pia ∼ 100 m) and externally by the dura membrane.At the center of the spinal cord lies the central canal, which in healthy subjects has a crosssectional diameter of approximately 0.5 − 2 mm.In the presence of hydro-/syringomyelia its diameter increases to approximately 3 − 4 mm.CSF is a water-like fluid with a density  CSF ≃ 1000 kg/m 3 and a kinematic viscosity  CSF ≃ 0.7 × 10 −6 m 2 /s.Its motion in the subarachnoid space is driven mainly by the intracranial pressure fluctuations that occur with each heartbeat as a result of the cyclic variation of the cerebrovascular blood volume (with frequencies on the order of  * CSF /(2) ∼ 1 Hz), resulting in a volume of CSF of approximately 1 ml that is cyclically pushed into and out of the spinal canal, as measured in multiple studies using magnetic resonance imaging [see, for example Coenen et al., 2019, Sincomb et al., 2022].The back-and-forth motion of CSF has an typical stroke length  * CSF ∼ 5 − 10 mm, and a characteristic velocity  * CSF ∼ 1 − 2 cm/s.The corresponding instantaneous Reynolds number is therefore on the order of  * CSF  * SSAS / CSF ∼ 50, and the corresponding Womersley number is on the order of  *

√︃
* CSF / CSF ∼ 8.Note that in the model problem to be presented below we will base the Reynolds number on the stroke length, so that its corresponding order of magnitude is  * CSF  * CSF / CSF ∼ 150.When the central canal is connected to the fourth ventricle, the fluid inside it is exposed to the same longitudinal oscillatory pressure gradient driven by the intracranial pressure fluctuations.Velocity measurements inside the central canal of healthy subjects are however not available in the literature because the diameter of the canal is too small for current MRI techniques.Nevertheless, there are a few measurements of flow inside large syringomyelia cysts, showing oscillatory motion in sync with that of CSF in the subarachnoid space surrounding the spinal cord [for example Honey et al., 2017, Luzzi et al., 2021].Various micro-anatomical features, such as nerve roots, denticulate ligaments, and trabeculae connect the pia membrane around the spinal cord with the outer dura mater.Among these, the trabeculae (Figure 1d) are especially relevant for the present problem, since they constitute a network of pillar-like structures that maintain the spinal cord in its position inside the spinal canal and in this manner give stiffness to the otherwise relatively flexible spinal cord tissue.The diameters of the trabeculae are on the order of  * trab ∼ 10 − 20 m, and they are spaced approximately  * trab ∼ 500 m apart [Stockman, 2006, Gupta et al., 2009, Salerno et al., 2020].In the problem formulation ( §2.1), we will show with an order of magnitude estimation of the governing fluid-structure interaction equation that the stiffness provided by the trabeculae is indeed the dominating effect that balances the normal force exerted by the two fluid layers on the separating spinal cord.
Inspired by the complex biofluid mechanical problem described in the previous paragraph, in the present paper we study a simplified, canonical configuration, consisting of two infinitely long layers of fluid that undergo an oscillatory motion parallel to a thin planar flexible wall that separates them.We investigate whether the separating wall will remain undeformed, the two oscillatory Stokes layers that form on both sides remaining unperturbed, or whether the configuration acquires a different oscillatory state in which the wall moves together with the two layers.We aim to answer this question by means of a Floquet linear stability analysis, exploring the influence of the different governing parameters such as the layer widths, the Reynolds number, and the wall stiffness.The study of fluid-structure interactions in simplified two-layer oscillating flow configurations such as the one presented here may prove helpful to elucidate some of the mechanisms that intervene in the complex, real-life flow described before.
The literature on the stability of time-periodic flows in the presence of compliant walls is limited.On the contrary, its rigid counterpart has been studied extensively Davis [1976].The stability of the planar Stokes layer near a rigid wall was investigated by Hall [1978] and Blennerhassett and Bassom [2002].Time-periodic channel and pipe flows with rigid walls have been studied by Yang and Yih [1977], von Kerczek [1982], Thomas et al. [2011], Pier and Schmid [2017], and Pier and Schmid [2021], using linear Floquet analyses, nonmodal stability analyses, as well fully nonlinear direct numerical simulations.Nevertheless, the absence of flexible walls limits the implications of their findings for the problem under consideration here.The hydrodynamic stability of flows near flexible walls has been mainly focused on steady configurations.Canonical examples include the Couette flow past a flexible membrane Kumaran and Srivatsan [1998], Thaokar and Kumaran [2002], Kumaran [2021], and the Poiseuille flow between compliant walls Carpenter and Garrad [1985], Davies and Carpenter [1997], Lebbal et al. [2022a].The coupling between the wall displacement and fluid velocity provides additional destabilizing mechanisms that enrich the stability characteristics drastically with respect to the rigid analogues of these flows.To the best of our knowledge, there are only a few studies in the literature on the stability of flows that contain both time-periodicity and flexible walls, among which the most relevant for the problem at hand work are the stability analyses of plane pulsatile channel flow between compliant walls by Tsigklifis and Lucey [2017] and Lebbal et al. [2022b], and the study of oscillating flow in a channel between a flexible wall and a rigid wall undergoing an in-plane pulsating motion by Thaokar and Kumaran [2004].Nevertheless, the main focus of these works is on the stabilizing or destabilizing influence of the modulation of the high-Reynolds-number mean flow, in contrast to the present work, in which the time-average of the basic oscillatory motion is zero.
The remainder of the paper is organized as follows.In §2, we describe the problem configuration, derive the basic oscillatory flow, and apply the Floquet formalism to the Navier-Stokes equations in combination with a simple model equation for the fluid-wall interactions to derive a set of stability equations that is to be solved numerically.The computational method employed to achieve this is set out in §3.The results are described in detail in §4.In particular, we focus on the marginal conditions for the onset of instability, investigating the influence of the perturbation wavelength and the control parameters of the problem, which are the layer widths, the Reynolds number, and the wall stiffness.Finally, concluding remarks are given in §5.

Governing equations and scaling
Based on the description of cerebrospinal fluid flow in the cervical spine given in §1, we consider here a simplified model problem that consists of an infinitely long channel of fluid with constant density  and kinematic viscosity , divided longitudinally into two layers of widths  * 1 and  * 2 by an infinitesimally thin flexible wall, as depicted in figure 2. Both layers are subject to the same oscillatory pressure gradient in the longitudinal  * -direction, −  * / * = Π * ℓ cos( *  * ) where  * denotes the angular frequency, which drives an oscillatory flow parallel to the walls.In the present analysis we are concerned with the stability of the basic oscillatory state in which the separating wall remains undeformed, and the motion in the two layers reduces to the unidirectional planar Womersley flow (see section 2.2).In the unstable regime, the separating wall undergoes deformations that are coupled to the fluid flow.
The choice of the particular fluid-structure interaction model considered here, i.e. a massless undamped spring-backed plate, is motivated by the cerebrospinal fluid flow on both sides of the spinal cord, the springs mimicking the trabeculae that keep the spinal cord in place within the spinal canal (figure 1d).The general form of the spring-backed plate-membrane equation, under the assumption of negligible tangential forces on the wall, and constraining the wall motion to the transverse -direction, takes on the form [see for example Carpenter and Garrad, 1985, Davies and Carpenter, 1997, Shankar and Kumaran, 2002] where ℎ * ( * ,  * ) is the transverse plate displacement, and [  * ] 1 2 is the pressure difference between the two sides of the membrane in the  * -direction.In the model,  * is the membrane mass per unit of area,  * is the wall-damping coefficient,  * is the flexural ridigity,  * is the longitudinal tension per unit width, and K * is the spring stiffness associated with the trabeculae.Taking into account that the spring stiffness of a single trabecular structure is on the order of  trab  * 2 trab / * SSAS , and that there are on the order of 1/ * trab trabeculae per unit area, the spring stiffness coefficient K * in equation ( 1) can be estimated to be Here  trab ∼ 12 kPa is the trabecular Young's modulus [Jin et al., 2011],  * trab ∼ 10 − 20 m the typical diameter of the trabeculae,  * SSAS ∼ 1 − 4 mm the width of the spinal canal (in other words, the length of the trabeculae), and  * trab ∼ 500 m the mean distance between trabeculae (figure 1d).Let us now compare the order of magnitude of the other terms on the left-hand-side of equation (1) to that of the spring-stiffness term, K * ℎ * .In these estimations, in addition to the aforementioned properties of the trabeculae, we need the Young's modulus and Poisson ratio of the pia mater,  pia ∼ 15 kPa and   ∼ 0.5, and its characteristic thickness,  * pia ∼ 100 m.We also assume small deformations in the transverse direction, on the order of ℎ *  ∼  * PIA , that span a longitudinal distance comparable to the stroke length of the cervical CSF motion,  * CSF ∼ 1 cm Carpenter and Garrad [1985], Howell et al. [2008].The relative orders of magnitude of inertia, flexural rigidity and longitudinal tension with respect to the string stiffness are, respectively, Damping is also considered to be negligible compared to the spring stiffness, with  *  * /K * ∼ 10 −2 Fiford and Bilston [2005].1It follows that in the fluid-structure interaction model (1), the spring stiffness is the dominating term on the left-hand-side that needs to balance the pressure difference between both sides of the wall induced by the fluid motion.Consequently, we focus our attention on the simplified model with  * =  * =  * =  * = 0, given by Note that this choice of fluid-wall coupling has been adopted in many biofluid mechanical problems in the literature Čanić et al. [2006], Westerhof et al. [2009], Sánchez et al. [2018].
When considering the fluid motion, it can be anticipated that in the regime of interest, i.e. near the critical conditions for the onset of instability, the thickness of the Stokes shear-wave layers near the walls is of the same order of magnitude or smaller than the channel half-widths, √︁ / * ≲  * 1,2 , while the characteristic oscillatory velocity in the bulk of the fluid layers is  * ∼ Π * ℓ /( * ), given by the balance between the driving pressure gradient and the local acceleration.Correspondingly, the typical length scale is the stroke length  * ∼  *  * −1 .Scaling lengths with  * , time with  * −1 , velocities with  * , and pressures with Π * ℓ 2 /( * 2 ), the dimensionless governing equations become (with the subindex  = 1, 2 distinguishing between layer 1,  < 0, and layer 2,  ≥ 0) The Reynolds number Re =  * 2 /( * ) =  *  * / = Π * ℓ 2 /( 2  * 3 ) and the dimensionless wall stiffness K = K * /Π * ℓ = K * /( *  * ) are the parameters that, together with the dimensionless channel widths  1 =  * 1 / * and  2 =  * 2 / * , govern the problem.Note that for the flow of CSF near the spinal cord suspended by the trabecular network inside the subarachnoid space, K ∼ 2 − 4, as can be obtained by introducing the values of the geometric and elastic parameters of the problem in 2. In §4 we will see that the fluid-structure interaction is indeed found to be optimal for order unity values of K.
The boundary conditions to be satisfied are the non-slip conditions at the channel walls and the kinematic condition at the deformable separating boundary, 1The reader is referred to §4.3 where we briefly explore the effect of the addition of this term, which only introduces little variations in the Floquet growth rate || and the critical Reynolds number Re  without changing the most unstable wavelength.Note that when  1 ,  2 become very large compared to the Stokes layer thickness, of order 1/

√
Re, the velocity profiles outside the Stokes layers tend to be uniform (sin , 0).In those cases, the Stokes layer on the nondeformable wall is so far away that it does not influence the flow near the flexible wall any longer.
In the present work we focus on Reynolds numbers of order unity or larger, so that the thickness of the Stokes layer is order unity or smaller.Therefore, a pragmatic approach to study the case  1 ,  2 → ∞ is to substitute the boundary conditions (11) by

Oscillatory base flow
The basic state about which the Floquet stability analysis is performed is the strictly parallel Womersleylike flow driven by the oscillatory pressure gradient when the separating boundary remains undeformed, given by ū2 with  = (1 + i) √︁ Re/2.Notice that, since the tangential stresses are not taken into account, the base flow of both streams must be in phase for the separating surface to remain undeformed.In the limit As an example, figure 3 shows the base flow for three different combinations of  1 ,  2 and Re.

Floquet linear stability analysis
Floquet theory Iooss and Joseph [2012] is employed to analyze the stability of the basic oscillatory state.
To that aim, small perturbations in the form of space-periodic Floquet modes are superposed onto the basic oscillatory state as ℎ = h() + ĥ() e i  e  .
(,  = ℎ, ) = û( = ℎ, ) + ĥ(, ) ū′  The eigenvalues  or, equivalently, the corresponding Floquet multipliers  = e 2   , determine the linear stability of the two-layer system.When, for a certain combination of governing parameters (K, Re,  1 ,  2 ) and for a certain wavenumber , all || < 1, the system is linearly stable to all perturbations of wavelength 2/.In the present work, we are concerned with obtaining the conditions for marginal stability, i.e. the values of (K, Re,  1 ,  2 ) and  for which a single value of  crosses the unit circle || = 1.

Numerical method
The system of stability equations ( 30)-( 34) is solved numerically using a finite number  of Fourier modes.A Chebyshev collocation method is used for the spatial discretization in the wall-normal (−)direction, employing a number  1 of collocation points in layer 1, and  2 points in layer 2. Since the -th Fourier mode is only coupled with modes  − 1 and  + 1, the general structure of the discretized eigenvalue problem is block-tridiagonal, where each set of square submatrices L ±1  and B ±1  , of size 3( 1 +  2 ) + 1, contains the discretized stability equations ( 30)-( 34).Their rows are arranged such that the first 3 1 rows encode the continuity and momentum equations ( 30)-( 32) for the fluid in layer 1, the next 3 2 rows encode the corresponding equations for the fluid in layer 2, and the last row contains the fluid-solid interaction equation (33).The boundary conditions (34) are implemented by replacing the corresponding rows in L ±1  and B ±1  .The numerical implementation was performed in MATLAB ® .The number of Fourier modes  and the number of discretization points  1 ,  2 was increased until numerical convergence was obtained (a Figure 5: The curve of marginal stability corresponding to || = 1 in the (Re, )-plane, for K = 1 and  2 =  1 = ∞.The critical Reynolds number Re cr , indicated with a dot, is defined as the minimum Re-value of the marginal curve.Note that the value Re cr = 25.9 corresponds to the critical condition shown in figure 4.
varying number of points was needed for varying layer depths  1 and  2 and varying Reynolds numbers Re).For all cases,  = 10 Fourier modes was found to be sufficient, whereas  1 and  2 ranged from 48 to 128.

Floquet eigenvalues and eigenmode analysis
The numerical solution of the eigenvalue problem that determines the linear Floquet stability of the two-layer configuration for a particular combination of canal widths  1 and  2 , stiffness K, Reynolds number Re, and perturbation wave number , yields a spectrum of eigenvalues , or, correspondingly, a spectrum of Floquet multipliers  = e 2   .The insets of Figure 4 show such spectra for K = 1,  1 =  2 = ∞, Re = 25.9 and three values of the wave number  = (0.8, 1.5, 2.2).Whereas for the smallest and largest wave numbers the most unstable mode has a nonzero imaginary part (left and right inset), there is a band of intermediate wave numbers for which the most unstable mode is purely real (middle inset), corresponding to perturbations of the flow field that oscillate synchronous with the base flow.It is this branch of Floquet modes that gets most amplified when the Reynolds number increases (compare in figure 4 the curve () for Re = 20 with that for Re = 25.9),crossing the unit circle || = 1 for  = 1.5 and Re = 25.9.A curve of marginal stability can be traced in the (Re − )-plane by obtaining, for each value of , the value of Re for which the most unstable mode has || = 1 (figure 5).We can then define a critical Reynolds number, Re cr , indicated in figure 5 with a dot, as the minimum Re-value of the marginal curve.Note that the value Re cr = 25.9 of figure 5 corresponds of course to the critical condition shown in figure 4. The condition (Re = 25.9,K = 1,  1 = ∞,  2 = ∞) thus lies on the marginal surface delimiting the linearly stable and unstable regimes in the parameter space spanned by Re, K,  1 and  2 .In section 4.2 different slices in this space are explored to see the individual influence of the different physical parameters on the stability of the system.
Figure 6 shows the time evolution of the perturbation eigenvector (a), together with that of the total, perturbed, flow (b), computed as the superposition of the perturbation eigenvector and the base flow, ( ū1,2 , v1,2 ) + 0.1( û1,2 , v1,2 )e i  e  , over the course of one oscillation cycle (with dimensionless period  = 2), for the marginally unstable case  = 1.52,Re = 26,  = 1,  1 = 3,  2 = 4.Note that the , where here † indicates the complex conjugate.Vectors indicate the direction of the flow and color represents the magnitude of the velocity perturbation in figure 6 (a), and of its superposition with the base flow in figure 6 (b).The perturbation corrugates the separating wall in a time-periodic fashion, with an associated spatiallyperiodic perturbed flow field that decays in the transverse direction  towards the upper and lower bounding walls.The spatial wave length is dictated by the wave number  as 2/ = 4.1, being  = 1.52 in this case.

Influence of the canal widths 𝐻 1 , 𝐻 2 and the wall stiffness K on the onset of instability
A parametric analysis has been carried out to study the influence of the widths  1 and  2 of the two fluid layers and the wall stiffness K on the linear Floquet stability of the system.Curves of marginal stability in the (Re, )-plane for different values of  1 are given for two illustrative cases:  2 → ∞, K = 1 (figure 7 a), and  2 = 1, K = 1 (figure 7 b).Reducing the canal width has a stabilizing effect, increasing the Reynolds number associated with the onset of instability over the entire range of wave numbers studied here.The most unstable wave number, corresponding to the critical Reynolds numbers indicated by the dots in figure 7, is not influenced appreciably by the canal widths, remaining close to  = 1.5 for the value K = 1 studied in figure 7, although a slight decrease of the most unstable wavenumber can be inferred in figure 7 (b) as  1 increases for  2 = 1.
Figure 8 slices the parameter space (Re,  1 ,  2 , K) by showing how the critical Reynolds number varies with  1 , for various values of  2 , and for four values of K, corresponding to the four panels (a)-(d).The stabilizing influence of decreasing the canal widths  1 and  2 becomes more pronounced the thinner the canals become.This is indicated both by the slope of the marginal curves near  1 = 1 and by the distance between the different marginal curves for different  2 .At this point, it is worth pointing out that, because of the symmetry of the configuration,  1 and  2 are completely interchangeable in the analysis.The stabilizing effect of confinement originates from the imposed attenuation of the velocity perturbations at the lateral walls (boundary condition 34).When unconfined ( 1 → ∞ and/or  2 → ∞), the eigenfunction decays exponentially in the transverse direction over a distance that scales with the perturbation wave length 2/ (stemming from the D 2 −  2 operator in the stability equations 31-32).This is illustrated in figure 9, which shows the velocity perturbation û2 (, ) at 8 instants of time over the course of one cycle for two critical cases, one case with a large wave number ( = 1.94,K = 0.5, Re = 30.1,figure 9 a), and another one with a small wave number ( = 0.76, K = 2.7, Re = 21.31,figure 9 b).As commented previously, when changing  1 or  2 , the most unstable wave number, corresponding to Re cr , hardly changes.Therefore, for smaller wave numbers (larger wave lengths), the transverse extent of the eigenfunctions being larger leads to the effect that decreasing  1 and/or  2 on the critical Reynolds number is noticed earlier (larger  1 or  2 ).On the contrary, for larger wave numbers (smaller wave   lengths), the eigenfunctions extend over a smaller distance, and the confinement is noticed later (smaller  1 or  2 ).This also explains why the influence of  1 and  2 occurs at different rates for the four panels of figure 8.Each panel corresponds to a different value of K and has a different associated wave number .
The occurrence of the instability and its associated wave number strongly depend on the value of the dimensionless stiffness K. Figure 10 shows, in blue, how the most unstable wave number  uns decreases with increasing stiffness K for symmetric channel configurations ( 1 =  2 ).Note that, in agreement with the results presented before, in the range 1 ≤ 1 =  2 < ∞ studied here,  uns (K) is essentially independent of  1 and  2 , and can therefore be presented in figure 10 as a single curve.The family of green curves in this figure represents the dependence of the Reynolds number, Re cr , with K. Except the case  1 =  2 = 1, there exist a value of K for which the Re cr is minimum.This minimum depends on the channel configuration, with the unbounded case being the most unstable, in the sense of exhibiting the lowest value of Re cr (K ≃ 3).The dimensionless wall stiffness K can be interpreted as the inverse of the square root of a reduced velocity, i.e.
√︁ 1/K = √︁  * /K * /( * −1 ), which is the ratio of the characteristic time associated with the spring-backed wall moving the adjacent fluid layer in the transverse direction, and the longitudinal oscillation time of the basic fluid motion.Therefore, optimal fluid-structure interaction can be expected for order-unity values of K, in agreement with the non-monotonic behavior of the marginal stability curves of figure 10.
In the limit of vanishing wall flexibility, i.e. when the spring stiffness K → ∞, the problem configuration reduces to that of oscillatory flow in a channel with rigid walls, which, for  1 ,  2 → ∞ further reduces to the oscillatory Stokes-layer flow over a rigid wall.The linear stability of the Stokes layer was first studied by Hall [1978], and revisited by Blennerhassett and Bassom [2002], who determined that the onset of instability occurs when the Reynolds number, based on the shear-layer thickness √︁ 2/ * , BB ≃ 10 6 and a wavelength  =  * / * =  BB /Re BB ≃ 2/270.Thus, the conditions for the wall deformation is optimally coupled to the oscillatory fluid motion, found here for K of order unity and Re ∼ 20 − 30 lie very far away from the critical conditions for the onset instability in the Stokes layer near a rigid wall.Capturing the transition from the present problem to its rigid counterpart requires a different treatment, using the shear-layer thickness instead of the stroke length as characteristic length scale, and lies out of the scope of the present work.

Influence of additional damping
As detailed and discussed in §2.1, the choice of the particular fluid-structure interaction model considered here, i.e. the massless undamped spring-backed plate, is motivated by the cerebrospinal fluid flow on both sides of the spinal cord in subjects with syringomyelia, the springs mimicking micro-anatomical features that keep the spinal cord in place within the spinal canal.Now, we briefly explore the effect of including a damping term in Eq. ( 10), where  =  *  * /Π * ℓ .This term is related with the amount of energy irreversibly removed from the wall.In the particular problem under consideration here, the order of the damping coefficient is estimated to be around  ∼ 10 −2 Fiford and Bilston [2005].Figure 11 shows the stabilizing effect of damping on the growth rate || of the perturbation (conditions similar to those of figure 4).Note that, the additional damping only introduces little variations in the Floquet growth rate || and the critical Reynolds number Re  , the most unstable wavelength remaining unchanged with respect to the undamped case.

Concluding remarks
We have investigated the stability of the fluid-structure interaction problem constituted by two layers of a fluid undergoing an oscillatory motion parallel to the compliant wall that separates them.The study of this canonical configuration is inspired by the cerebrospinal fluid flow occurring in hydro-/syringomyelia, a pathology of the central nervous system characterized by the accumulation of fluid within the spinal cord, and its aim is to shed light on the underlying mechanical processes.The problem is governed by the Navier-Stokes equations for a Newtonian, incompressible, fluid, together with a constitutive equation for the compliant solid, for which a spring-backed plate equation was used, only accounting for the stiffness K and assuming wall-normal deformations.We investigated the linear stability of the problem by means of a Floquet analysis, focusing on the effect of the different control parameters of the problem, namely the dimensionless fluid layer widths,  1 and  2 , the Reynolds number, Re, and the dimensionless wall stiffness, K.
Around criticality, i.e. when the absolute value of the Floquet multiplier of the most unstable mode crosses the unit circle, it was found that the perturbations of the flow field oscillate synchronous with the base flow.The associated dimensionless perturbation wavenumber  is of order unity, in other words, the spatially periodic deformations of the separating wall have a length 2/ * that is of the order of the stroke length  * of the oscillatory fluid motion.The exact value of that wavenumber depends strongly on the wall stiffness K, as does the critical Reynolds number Re cr , underscoring the role that K plays as a reduced velocity, comparing the characteristic time of the spring-induced transverse motion with the longitudinal fluid oscillation period.Our results also reveal that reducing the channel widths has a stabilizing effect, i.e. a larger Reynolds number is necessary to trigger the onset of instability over the entire range of perturbation wave numbers herein considered.In particular, such stabilizing effect has been found to become stronger as the fluid layers become thinner.The latter stems from the attenuation of the velocity perturbation close to the lateral rigid walls.
The present work aims to be a first step in understanding the complex fluid-structure interaction problem encountered in syringomyelia.It is encouraging that we find the most unstable wave lengths to be approximately 2/1.5 ≃ 4 times the stroke length  * , which, with  * CSF ≃ 1 cm in the spinal canal, yields 4 cm, comparable to the size of fluid accumulation.This may suggest that syrinx sizes lie in the range in which optimal fluid-structure interaction is achieved, and in which the motion of the intrasyringal fluid is most easily coupled with that of the cerebrospinal fluid surrounding the spinal cord.Future efforts will be focused on the validation of the results using laboratory experiments and direct numerical simulations.

Figure 1 :
Figure 1: (a) Schematic overview of the main features of the anatomy of the cervical spine.b) Schematic representation of syringomyelia associated with Type I Chiari Malformation.c) Schematic representation of hydromyelia.d) Diagram and photograph Cloyd and Low [1974] of the micro-anatomy inside the spinal subarachnoid space; here AM: arachnoid mater, AT: arachnoid trabeculae, A: artery, SAS: subarachnoid space, DM: dura mater, PM: pia Mater, SPC: spinal cord.

Figure 2 :
Figure 2: Schematic representation of the flow configuration.

Figure 4 :
Figure 4: Floquet growth rate || as a function of the perturbation wave number , for two values of the Reynolds number, Re = 20 and Re = 25.9, and K = 1,  1 =  2 → ∞.Insets show the spectra of Floquet multipliers  for three different wave numbers.The critical conditions correspond to Re = 25.9 and  = 1.5.

Figure 7 :
Figure 7: Curves of marginal stability in the (Re, )-plane, for various values of  1 for the cases ( 2 = ∞, K = 1) (a) and ( 2 = 1, K = 1) (b).Solid dots indicate the critical Reynolds number Re cr as the minimum Re-value of each marginal curve.

Figure 8 :
Figure 8: Variation of the critical Reynolds number Re cr with  1 , for different values of  2 , and four values of K: (a) K = 0.5, (b) K = 1, (c) K = 2 and (d) K = 2.7.Note that the wave number is approximately  ≃ 1.94 for all cases in (a),  ≃ 1.5 for all cases in (b),  ≃ 0.98 for all cases in (c), and  ≃ 0.76 for all cases in (d).

Figure 11 :
Figure 11: Floquet growth rate || as a function of the perturbation wave number , for Re = 26, and K = 1,  1 =  2 = ∞, when, on top of the spring stiffness, a damping term is included in the fluid-wall interaction model, for three values of the corresponding dimensionless damping coefficient,  = 0,  = 0.02,  = 0.04.