Phase-reduction for synchronization of oscillating flow by perturbation on surrounding structure

Abstract Regulation of fluid flow by deformations of the surrounding elastic structure is observed in many natural and artificial system, such as in the cardiovascular system. As the first step to study the regulation of oscillating flows, we consider synchronization of vortex shedding past a cylinder within an elastic structure with a sinusoidal external forcing. We use phase-reduction theory to evaluate the synchronization characteristics of the oscillating fluid–structure coupled dynamics. We find that the phase-sensitivity function, which characterizes the phase-response of the oscillation, is significantly affected by the Cauchy number and slightly affected by the fluid-to-structure density ratio and Poisson's ratio of the structure material, for fixed model configuration and Reynolds number. The predicted synchronization characteristics are in close agreement with results from direct numerical simulations. The synchronization region is maximized when the sinusoidal perturbation is applied near the downstream end of the cylinder. These findings open further possibility for the utilization of phase-reduction theory to characterize synchronization in other practical problems exhibiting fluid–structure coupled dynamics, such as in biological systems and the control of microfluidics.

cardiac cells (Kojima, Kaneko & Yasuda 2006) and vortex shedding (Zdravkovich 1996). Such systems can be interpreted as limit-cycle oscillators. One interesting phenomenon observed in oscillating dynamics is synchronization (Kuramoto 1984;Pikovsky, Rosenblum & Kurths 2001;Ermentrout & Terman 2010). Although such dynamical systems are generally nonlinear and can be infinite-dimensional, their periodic nature allows the dynamics to be represented in a simplified manner using phase-reduction theory, which is a powerful tool that allows synchronization characteristics to be analysed by considering only the phase dynamics of the system (Kuramoto 1984;Pikovsky et al. 2001;Ermentrout & Terman 2010;Nakao 2016).
Phase-reduction theory has been extended for many different kinds of oscillatory dynamics such as in time-delayed systems and hybrid systems (Kotani et al. 2012;Novičenko & Pyragas 2012;Shirasaka, Kurebayashi & Nakao 2017). In the field of fluid mechanics, recent work by Kawamura & Nakao (2013, 2015 considers the extension of phase-reduction theory for periodic spatiotemporal patterns to analyse Hele-Shaw flows. More recently, a computational study for the lock-on of vortex shedding to oscillatory actuation was considered in the work of Taira & Nakao (2018) and Khodkar & Taira (2020).
In cases where the oscillating fluid flow comes into contact with an elastic structure, the dynamics need to be considered as a problem of fluid-structure interaction (FSI). The surrounding structure displacement will cause movement of the fluid boundaries. Conversely, the fluid flow exerts force on the fluid-structure interface, causing structural displacements (Richter 2017). A common example of an oscillatory phenomenon involving FSI is a vibrating flexible structure within vortex shedding (Gomes & Lienhart 2013). Self-oscillation induced by FSI is also observed in the study of micro-electromechanical systems (Ducloux et al. 2007). In biological systems, it is known that the efferent activity of the sympathetic nervous system on vascular smooth muscle alters the contractility of blood vessels, which in turn modulates the periodic blood flow (Kotani et al. 2005;Shiogai, Stefanovska & McClintock 2010). Theoretical formulation using phase-reduction for elastohydrodynamic synchronization of beating flagella has been considered in the work of Kawamura & Tsubaki (2018). However, much remains unknown about oscillation and synchronization in fluid-structure coupled dynamics.
In this work, we present the utilization of phase-reduction for FSI dynamics as the first step to understand how oscillatory flow can be regulated by weak perturbation on the surrounding structure. We describe the model problem of vortex shedding past a cylinder within an elastic structure as well as the phase-reduction theory in § 2. The dependence of the phase-response property on the material types and the location of the perturbation as well as the synchronization characteristics are discussed in § 3. Finally, concluding remarks are given in § 4. The method presented here can be seen as a simplified yet powerful approach to analyse a complex dynamical system involving multiphysics, such as FSI, by reducing it into single scalar phase dynamics.

Model definition
In this study, we define a model of two-dimensional laminar incompressible flow past a cylinder within a channel surrounded by an elastic structure, as shown in figure 1. The cylinder itself is considered as a static rigid body with diameter d = 0.1 m centred at (x, y) = (0.2 m, 0.1 m). The fluid and structure are fully coupled (Horn & Turek 2006; x (m) Richter 2017). Numerical simulations are conducted based on the finite element method using COMSOL Multiphysics 5.5. The fluid flow is governed by the incompressible Navier-Stokes equations in the arbitrary Lagrangian-Eulerian (ALE) formulation (Quarteroni, Tuveri & Veneziani 2000;Horn & Turek 2006) such that where ρ f , μ, p, u f , u m and I represent fluid density, dynamic viscosity, pressure, fluid flow velocity, spatial coordinate velocity and identity matrix. The external volume forces are assumed to be zero. The fluid flows inside the channel from left to right. The outlet is prescribed with zero-pressure boundary condition. A no-slip boundary condition is prescribed at the fluid-structure interface and around the cylinder wall. A parabolic velocity profile is prescribed at the inlet such that u(0, y) =Ū f 6y(H − y)/H 2 , whereŪ f is the mean inlet velocity and H = 0.2 m is the channel width. The surrounding structure is modelled as an isotropic linear elastic material. The dynamics is governed by where ρ s , w s , α, β and P represent structure density, displacement, mass damping coefficient, stiffness damping coefficient and the first Piola-Kirchoff stress tensor (Formato et al. 2019). The stress tensor P is characterized by Young's modulus E and Poisson's ratio ν of the structure material according to Hooke's law. In this work, we define α = 0.2 s −1 and β = 0.1 s. The width of the elastic structure is set at 0.02 m for both the top and bottom sides. Zero displacement is prescribed at the left and right edges of the structure. The perturbation εη s (x, t) is defined as a localized force at the upper boundary of the structure, as shown in figure 1, such that where x is the spatial coordinate, δ(x − x 0 ) describes the spatial impulse profile andê y is the unit vector in the y direction, while the negative sign signifies that the perturbation is applied in the direction into the structure domain. The time-varying component η(t) is defined as sin(ω f t), for sinusoidal perturbation. (2.4) Both the spatial and temporal impulse profiles are approximated as Gaussian functions with σ x = 0.01 m and σ t = 0.01 s. We use the impulse perturbation to evaluate the phase-sensitivity function of the oscillating flow, which is necessary in the phase-reduction analysis, and the sinusoidal perturbation to evaluate synchronization of the oscillating flow (more later in § 2.2). The perturbation amplitudes are chosen such that ε = 0.001 N for the case of impulse perturbation and ε ∈ [0, 25] N for the case of periodic perturbation. The kinematic and dynamic coupling conditions are satisfied at the fluid-structure interface Ω (Richter 2017), that is, where u s is the velocity of the structure displacement, F n is the normal force per unit area exerted on the structure by the fluid andê n is a unit vector normal to the fluid-structure interface. We characterize the model problem by the cylinder-diameter-based Reynolds number Re, Cauchy number C Y and fluid-to-structure density ratio M (DeNayer et al. 2018), which are defined as (2.6a-c) Table 1 shows the material properties that are considered in this study. Using properties of material type 1, but by modifying the value of μ, it is found that limit-cycle oscillation occurs for 90 ≤ Re ≤ 300. For all material types defined, the fluid properties are chosen such that Re = 100. Material type 1 is used as base values for comparisons. For types 2 and 3, the value of C Y is modified by changing E with a ratio of 1.25 times the base value. For types 4 and 5, the value of M is modified by changing ρ s with a ratio of 1.25 times the base value. For types 6 and 7, ν is modified with a ratio of 1.2 times the base value to keep it within the range for an isotropic linear elastic material (−1.0 < ν < 0.5). For type 8, both fluid and structure properties are changed while keeping the values of Re, C Y , M and ν equal to the base values.
The computational domain is discretized using a distributed quadrilateral mesh with a total number of 2864 elements. The elements around the cylinder and the fluid-structure interface are more refined with maximum size limit at 6.72 × 10 −3 m and minimum size limit at 9.6 × 10 −5 m. Implicit time stepping is used and the simulation output is stored every t = 0.0005 s.

Phase-reduction
Here, we briefly describe the phase-reduction theory and its use for analysing synchronization characteristics (Nakao 2016;Taira & Nakao 2018). Let us write the governing dynamics of an oscillating FSI dynamics that is perturbed by an external force as where X = [u f , w s ] is the system state, F {X } is the system dynamics and the perturbation εη(x, t) is sufficiently small. We assume that the system has an exponentially stable limit-cycle solution X 0 with frequency ω n , such that X 0 ( that maps the system state to a scalar phase value θ ∈ [0, 2π) can be introduced (Nakao, Yanagita & Kawamura 2014), such that the phase of the system θ = Θ[X ] always increases with a constant frequency aṡ θ(t) = (δΘ/δX ) · F {X } dx = ω n in the basin of the limit cycle when the perturbation is absent, where δΘ/δX is the functional derivative of Θ[X ] at X = X (x, t). At the lowest-order approximation (neglecting the higher-order-terms), δΘ/δX can be approximately evaluated at X = X 0 (x, θ/ω n ) on the limit cycle and the phase dynamics under weak perturbation can be found aṡ is known as the phasesensitivity function. If Z (x, θ) is known, the influence of any perturbation function εη(x, t) can be determined through (2.8). In particular, if the perturbation on the structure is spatially localized in the form η( In this work, we evaluate Z y (θ ) using the direct impulse perturbation approach (Ermentrout & Terman 2010;Nakao 2016), with the perturbation function εη s (x, t) defined in (2.3a,b) and (2.4). By introducing the impulse at certain phase θ (determined by t 0 ), the phase dynamics will be affected, where phase shift will be observed along the limit-cycle orbit. This asymptotic phase shift is known as the phase-response function g(θ ; εê y ), which is dependent on the amplitude, location, direction and phase where it is introduced. Hence, the phase-sensitivity function at phase θ can be evaluated as (2.9) Therefore, we can evaluate Z y (θ ) by introducing impulse perturbations over the entire range of θ. By evaluating g(θ ; εê y ) using appropriate observables showing periodic oscillations, Z y (θ ) can be determined using a limited number of measurements. In this study, we consider the oscillating lift coefficient C L on the cylinder boundaries as the observable. The relative phase value of θ = {0, 2π} is defined at the minimum value of C L oscillation. The representation in phase dynamics allows further analysis to determine the conditions in which the original dynamics will synchronize to a periodic perturbation with period T f and frequency ω f = 2π/T f . We consider the phase difference φ(t) between the phase of the original oscillatory dynamics θ(t) and the periodic perturbation ω f t such that (2.10) Combining with (2.8) for the case where a perturbation is applied on the structure and by applying the averaging approximation (Ermentrout & Terman 2010), we can evaluate the rate of change of φ(t) such thaṫ where Γ (φ) is known as the phase-coupling function, which is 2π-periodic, and ψ = ω f t. By conducting stability analysis on the phase dynamics in (2.11a), we can determine the criteria for which the original oscillating dynamics will synchronize to a periodic perturbation, that is, when |φ| → 0 and ω n converges to ω f . In this case, the original limit-cycle oscillation exhibits phase-locking to the periodic perturbation asymptotically at a stable fixed point φ. Through observation of (2.11a), we can find that the phase dynamics stability is satisfied if (2.12) Using this stability criterion, we can determine the boundaries of the synchronization region over the (ω f /ω n )-ε space. These synchronization boundaries are also known as the Arnold tongue (Pikovsky et al. 2001) and will be discussed later in § 3. On the contrary, periodic phase slip will occur if |φ| > 0 and the phase difference φ(t) will continuously increase over time (Nakao 2016). We can define the frequency of this periodic phase slip as f slip = 1/T slip , where T slip is the interval of the periodic phase slip, which is calculated as a function of frequency difference ω n − ω f . When the original oscillating dynamics is phase-locked to the periodic perturbation, the mean frequency of the oscillation is equal to ω f and f slip ≈ 0. Outside of the phase-locking region, phase slip will occur and f slip can be calculated as . (2.13a,b) Depending on the frequency difference, the approximated value of f slip can be positive or negative. The characteristics of the f slip curve over the span of ω n − ω f can be compared with the results from actual direct numerical simulations (DNS), as seen later in § 3.  π/2 π 3π/2 2π 0 π/2 π 3π/2 2π π/2 π 3π/2 2π

Phase-sensitivity analysis
In what follows, we compare the phase-sensitivity function Z y (θ ) for different material types, given that the impulse perturbation defined in (2.3a,b) and (2.4) is introduced. The findings in this section can be further utilized to evaluate the phase-coupling function and to analyse the synchronization characteristics. For each material type, 11 actual simulations were performed at equal phase intervals over θ (shown by markers in figures 2 and 3) and the in-between values are obtained by the modified Akima interpolation method using MATLAB. Figure 2(a) shows a comparison of the C L oscillations between material types 1 and 8 in which both the fluid and structure properties are changed but all the characteristic numbers are kept equal. The period of the C L oscillations is found at T n = 0.403 s for material types 1 to 7, while for material type 8 the period is found at T n = 0.323 s, showing that the oscillation is mainly characterized by the fluid properties (see table 1). For all material types, the vortex shedding frequency is 3.5 times or more higher than the dominant harmonic frequency of the structure, such that large structural resonance due to the vortex shedding does not occur. Note that in the present study damping terms are included in the structural dynamics, such that the structure harmonic modes have negligible effect on the limit-cycle oscillation as time goes to infinity.
We confirm that the chosen characteristic numbers used in this study determine the phase dynamics by the results shown in figure 2(b), where close agreement can be seen between Z y (θ ) evaluated for material types 1 and 8 (shown as solid and dotted lines, respectively). As shown in figure 2(a), the C L oscillation periods for material types 1 and 8 differ due to the different mean inlet velocities. Nevertheless, when all the values of Re, C Y , M and ν are matched (see table 1), the obtained phase-sensitivity functions are in close agreement. This shows that for an equal model configuration, the phase dynamics is characterized by the values of Re, C Y , M and ν. Figure 2(b) also shows that the phase-sensitivity function Z y (θ ) is significantly affected by the perturbation location. We compare the results for perturbations applied in between the upstream and downstream sides relative to the cylinder centre on the x axis at x = 0.2 m. Smaller mean and peak-to-peak values are observed for perturbations applied at the upstream side (x 0 < 0.2 m) while larger values are observed for the cases where perturbations are applied at the downstream side (x 0 > 0.2 m). In the upstream region, flow separation has yet to occur such that the oscillating flow is not fully developed. Hence, for an equal impulse perturbation amplitude, the phase-response is considerably smaller than in the case when the perturbation is applied at the downstream side where vortices are starting to develop. These also suggest that the required periodic perturbation amplitude and frequency to achieve synchronization, as defined in (2.12), are also dependent on the perturbation location, which will be studied in the next subsection.
We further investigate the effect of each characteristic number by comparing the response to impulse perturbation introduced at x 0 = 0.25 m. As mentioned in § 2.1, we modify the Cauchy number C Y and the fluid-to-structure density ratio M by changing the Young's modulus E and structure density ρ s , respectively. Figure 3(a) shows the comparison for material types 1, 2 and 3, where Z y (θ ) is compared for different values of C Y . The mean and peak-to-peak values of Z y (θ ) increase as the value of C Y increases. Figure 3(b) shows the comparison for material types 1, 4 and 5, where Z y (θ ) is compared for different values of M. The peak-to-peak value of Z y (θ ) increases as the value of M increases while the mean value is unaffected. Figure 3(c) shows the comparison for material types 1, 6 and 7, where Z y (θ ) is compared for different values of Poisson's ratio ν. The mean and peak-to-peak values of Z y (θ ) increase as the value of ν decreases.
Comparisons on the effect of the structure material and fluid-structure coupling characteristic numbers show that the changes in M and ν have smaller effects on the phase-sensitivity function than that of C Y , even taking into account the narrower range of ν due to the parameter restriction. The change in C Y significantly affects the mean value of the phase-sensitivity function rather than its waveform. Although the effects of C Y , M and ν are different, these values characterize the overall phase-response to the perturbation on the elastic structure.
Observation of figures 2 and 3 shows that the impulse perturbation causes phase advancement to the oscillating flow for all the chosen material types and perturbation locations. All the resulting phase-sensitivity functions have positive non-zero mean values. Similar phase-response is typically found in class I neurons (Ermentrout & Terman 2010). These suggest that, if we periodically push the upper boundary of the structure, then the perturbation frequency must be higher than the natural frequency of the oscillating flow in order to achieve synchronization. Such limitation can be avoided if we consider that the perturbation is in the form of a periodic pushing and pulling action on the upper boundary, which is defined in this work as a zero-mean sinusoidal function shown in (2.4).

Synchronization analysis
We consider material type 1 and the periodic sinusoidal perturbation defined in (2.3a,b) and (2.4) being used for the following synchronization analysis. Given that the phase-sensitivity function Z y (θ ) can be determined, we can evaluate the phase-coupling function Γ (φ) using (2.11b) for the given perturbation function. The phase-coupling functions evaluated at different perturbation locations x 0 are shown in figure 4(a). As can be deduced from the peak-to-peak values of Z y (θ ) shown in figure 2(b), the evaluated phase-coupling functions show largest variation over φ when the perturbation is applied at x 0 = 0.25 m. Variation of Γ (φ) over negative and positive values suggests that a perturbation frequency lower than the natural frequency of the oscillating flow can also lead to synchronization.
Several other numerical simulations for perturbation locations x 0 ∈ [0.1, 0.4] m are conducted to compare their synchronization regions, as shown in figure 4(b). The markers indicate locations where simulations are actually performed. By the criterion for synchronization in (2.12), larger variation of Γ (φ) will result in a wider synchronization region. For the given model problem, we found that the widest synchronization region is achieved at 0.30 < x 0 < 0.35 m. In general, a wider synchronization region is achieved given that the perturbation is applied near the downstream end of the cylinder in which the vortices are fully developed. The synchronization region shrinks as the perturbation location goes farther downstream from the cylinder as the vortices start to decay.
In what follows, we consider the case where periodic perturbation is applied at x 0 = 0.25 m. Using the synchronization criterion in (2.12), we can construct a representation of the synchronization region in the form of an Arnold tongue, as shown in figure 4(c). We compare the approximated synchronization region with results obtained from DNS. For the actual simulations, the C L oscillation and the periodic perturbation are considered to be synchronized when their phase difference φ stops growing and exhibits small fluctuations around a certain value (i.e. phase-locking occurs). It can be seen that the approximated synchronization region agrees well with the results from actual simulations of the nonlinear FSI dynamics, especially for smaller values of the perturbation amplitude ε. The effect of nonlinearities become more apparent for ε ≥ 15 N in which the actual dynamics starts to deviate from the approximation obtained using phase-reduction theory. As also noted in the work of Taira & Nakao (2018), the identification of the Arnold tongue through the phase-coupling function is attractive since we only need to conduct a small number of numerical simulations over the phase θ ∈ [0, 2π) to determine the synchronization region.
The elastic structure displacement oscillation in the y direction at x = 0.25 m has a mean value of −2.77 × 10 −6 m and a peak-to-peak value of 0.54 × 10 −6 m when the periodic perturbation is not applied. For the synchronized case where a periodic perturbation with amplitude ε = 25 N is applied, the displacement has a mean value of −7.97 × 10 −6 m and a peak-to-peak value of 0.004 m. It is also observed that the oscillating displacement in the y direction of the elastic structure at the bottom side synchronizes with the periodic perturbation due to the fluid-structure coupling.
We compare the approximated f slip characteristics evaluated using (2.13a,b) with the results obtained from DNS, as shown in figure 4(d). In this work, the f slip characteristic is compared for perturbation amplitude ε = 10 N. However, comparisons can be conducted for any value of ε. The approximated phase-locking region is shown between the dashed vertical lines. We can see that the f slip characteristics obtained from DNS (shown with the markers) is in close agreement with the approximated characteristics obtained from the calculations. However, the actual phase-locking region is slightly biased towards the lower perturbation frequency. This also agrees with the findings seen in figure 4(c) in which the actual synchronization region is biased towards the lower perturbation frequency, especially for larger perturbation amplitudes. These findings further confirm that the linear approximation using phase-reduction captures the synchronization characteristics of the original nonlinear oscillating dynamics.
To the best of our knowledge, this is the first study that uses phase-reduction theory to determine synchronization characteristics of an oscillating flow in FSI dynamics. Here, we have presented the idea of how synchronization between oscillatory flow and periodic forcing on the surrounding elastic structure can be achieved, as well as the use of phase-reduction theory to determine the region of synchronization.
A recent study in the field of microfluidics (Sun et al. 2017) presented a design of microfluidic oscillators where self-sustained flow oscillation induced by impinging jet flow on a bluff body is incorporated for mixing of fluids . It is also known that fluid mixing can be induced by vibration on the wall of the fluid container (Carlsson, Sen & Löfdahl 2005). Combination of these two concepts might reveal an alternative method to enhance the mixing performance, and phase-reduction theory can be applied to estimate the required parameters for obtaining the desired flow regulation, such as presented in this study.
The understanding of periodic fluid flow regulation by the motion of an elastic structure is also important for studies of blood pressure regulation due to the movement of vascular branches induced by the activity of the sympathetic nervous system (Kotani et al. 2005).
There are, of course, other synchronization cases that can be considered, such as periodic perturbation applied on a structure submerged inside an oscillating fluid flow. In addition, the approach using phase-reduction theory can also be extended to analyse the synchronization around harmonic frequencies. Recent experimental work and mathematical modelling (Barros et al. 2016;Herrmann et al. 2020) show how symmetric forcing leads to subharmonic synchronization while antisymmetric forcing leads to harmonic synchronization. These phenomena may also be observed in FSI dynamics as well, and phase-reduction theory can also be utilized to uncover synchronization characteristics under symmetric and antisymmetric forcing.

Conclusion
We have applied phase-reduction theory to find the synchronization characteristics of vortex shedding by periodic perturbation on the surrounding elastic structure. We have conducted parametric studies to observe how several dimensionless parameters known in FSI dynamics affect the phase dynamics of the limit-cycle oscillation. It is found that the phase-sensitivity function is significantly affected by the Cauchy number, whereas it is affected only rather slightly by the fluid-to-structure density ratio and Poisson's ratio, given that the model configuration and fluid flow characteristic defined by the Reynolds number are equal. The perturbation location on the elastic structure also significantly affects the phase-sensitivity function. Furthermore, it is confirmed that the estimated synchronization characteristics are in close agreement with the results obtained from DNS by comparing the synchronization region and the periodic phase slip characteristics. The synchronization region is maximized given that the sinusoidal perturbation is applied near the downstream end of the cylinder. Our findings open the further possibility for the utilization of phase-reduction theory for synchronization analysis in other practical problems involving oscillations in dynamics governed by fluid-structure interaction, such as in biological systems and control of microfluidics, using only a limited number of experiments or simulations. Declaration of interests. The authors report no conflict of interest.