Controlling fluidic oscillator flow dynamics by elastic structure vibration

In this study, we introduce a design of a feedback-type fluidic oscillator with elastic structures surrounding its feedback channel. By employing phase reduction theory, we extract the phase sensitivity function of the complex fluid–structure coupled system, which represents the system’s oscillatory characteristics. We show that the frequency of the oscillating flow inside the fluidic oscillator can be modulated by inducing synchronization with the weak periodic forcing from the elastic structure vibration. This design approach adds controllability to the fluidic oscillator, where conventionally, the intrinsic oscillatory characteristics of such device were highly determined by its geometry. The synchronization-induced control also changes the physical characteristics of the oscillatory fluid flow, which can be beneficial for practical applications, such as promoting better fluid mixing without changing the overall geometry of the device. Furthermore, by analyzing the phase sensitivity function, we demonstrate how the use of phase reduction theory gives good estimation of the synchronization condition with minimal number of experiments, allowing for a more efficient control design process. Finally, we show how an optimal control signal can be designed to reach the fastest time to synchronization.

Self-sustaining oscillation can be observed as a natural phenomenon in the study of fluid physics. Vortex shedding 1 is one prime example of such phenomena, and is well-known to potentially cause flow-induced vibration 2 . On the contrary, fluidic oscillators are devices designed to generate oscillating flow which can be exploited for other uses. Ever since its invention, fluidic oscillators have found many uses such as flow separation control, drag reduction, flowmeter design, and fluid mixing [3][4][5][6][7] .
One of the most attractive characteristics of the fluidic oscillator is that it requires no moving component. In the case of the feedback-type fluidic oscillator, the self-sustaining oscillation is generated because the fluid jet supplied into the device will be attached to one side of the inner walls due to Coanda effect, causing a small fraction of the fluid jet to return through the feedback channel and deflect the main inlet flow to the other side. This process repeats alternately, creating a self-sustaining oscillation [8][9][10] , as illustrated in Fig. 1. However, this characteristic imposes a limitation in which the geometry of the device defines the resulting characteristics of the oscillatory flow, including its frequency 10,11 . The simplest method to modulate the oscillation frequency is to change the inlet pressure, which will consequently increase the inlet flow rate. However, an increase in inlet pressure may not be desirable as the device must be able to withstand larger pressure value 12 . Furthermore, the change in inlet flow rate may also change the flow regime from laminar to turbulent flow.
Mechanisms to generate oscillatory flow in a fluidic oscillator with frequency independent of its inlet flow condition has been considered in past studies. Methods that employ direct excitation such as by a piezoelectric actuator 13 or a plasma discharge mechanism 14 are common examples. However, the fluid flow itself does not intrinsically exhibit oscillatory dynamics in such cases. Other methods use the generated sweeping jet as additional feedback 9 , or use a cascaded oscillator geometry 15,16 .
In order to analyze a fluidic oscillator performance and its resulting fluid flow characteristics, a complicated experimental work is usually required, from manufacturing the device to setting up the equipment necessary for data acquisition, such as with hot-wire anemometry 3,15,16 . Numerical simulation allows calculation of the fluid flow characteristics of a fluidic oscillator design before the actual experiment. However, apart from conducting multiple simulations with different parameters, it is generally difficult to estimate the controllability of the system or to find optimal settings under given restrictions.
On the other hand, the use of reduced order models to describe oscillatory flow has become popular as it allows for a simplified description of the underlying physics, and predicts the resulting fluid flow characteristic for www.nature.com/scientificreports/ a given control input [17][18][19][20][21] . For oscillating flow control, periodic forcing to induce synchronization is a rigorously studied subject [21][22][23] . Applying model reduction theory will benefit this field by simplifying its analysis. Model reduction based on phase reduction theory was first studied by Taira and Nakao, who found that low-dimensional dynamic representation using a single phase equation could estimate synchronization properties of a vortex shedding past a cylinder 24 . While the method is more commonly used in studies of biological systems 25,26 , it has gained much interest in the fluid physics community in the recent years 22,27,28 . Other studies have expanded its usage to fluid-structure interactions 29 and thermoacoustic systems 30 .
In this work, we propose a control mechanism to modulate oscillation frequency of a feedback-type fluidic oscillator by synchronization with weak external forcing in the form of vibration of elastic structures attached to the fluidic oscillator body. In our proposal, the intrinsic oscillatory flow is maintained even without external forcing, while the moving component adds controllability to the system. Phase reduction theory is used to simplify the analysis of the synchronization condition of the complex fluid-structure interaction dynamics, and allows estimation of the control region. We demonstrate how different placement of the elastic structure leads to a different phase response and analyze the change in fluid flow properties both for when synchronization is achieved or not. Our results show that the synchronization condition is well estimated using phase reduction theory, and that the synchronization induces changes in physical properties of the oscillating fluid motion. In addition, we demonstrate how an optimal control framework can be used to design the shape of the periodic forcing signal that allows for faster synchronization time.

Results and discussion
Oscillatory characteristics. We begin by describing the design of the fluidic oscillator model used in this study and observing its oscillatory characteristics. We use a typical feedback-type fluidic oscillator geometry in our simulation setup, as shown in Fig. 2. This geometry is adapted from a previous study by Seo et. al. 31 with the addition of an elastic structure around the fluidic oscillator being the most prominent difference. In our setup, a laminar incompressible fluid flow is considered and the elastic material is isotropic with linear elasticity. See the "Methods" section for a detailed explanation of the simulation setup. Two observables were measured to evaluate the oscillatory characteristics. The first one is the pressure loss coefficient C P , which correlates to the energy loss required to maintain the oscillation. A smaller value of C P corresponds to better energy efficiency. The second one is the driving force coefficient C D , which correlates to the force driving the flow to the feedback channels. Hence, a larger value of C D suggests a larger lateral momentum inside the mixing chamber of the fluidic oscillator. These two observables are measured as where p 1 -p 4 are the dimensionless pressure values measured at four specific locations in the computational model, d i is the power nozzle width, and d is the feedback channel width, as shown in Fig. 3. The time series for both observables are shown in Fig. 3a. It is clear that C P is oscillating twice as fast as C D because the sweeping fluid motion passes through the outlet twice in a single sweep. The time course of C D is obtained by measuring the pressure difference in the lateral direction, corresponding to the oscillatory characteristics of the lateral sweeping motion inside the fluidic oscillator, which is of our interest. Hence, the time course of C D will be used as the basis for the phase reduction analysis.
Since the addition of the elastic structure is intended to give controllability of the oscillatory dynamics, the intrinsic oscillatory dynamics inside the fluidic oscillator is expected to be similar to a conventional fluidic www.nature.com/scientificreports/ oscillator with no moving parts when there is no periodic forcing applied to the elastic structure. To confirm this, we compare the resulting Strouhal number St for a given Reynolds number Re with and without the addition of the elastic structures on the design. We confirm that the addition of the elastic structures itself does not affect the intrinsic oscillatory dynamics of the fluid flow inside the fluidic oscillator, as shown in Fig. 3b (see "Methods" section for details). Since the Strouhal number can be considered as the normalized oscillation frequency with respect to the fluid velocity, Fig. 3b shows that the oscillation frequency scales linearly with the Reynolds number, which is an expected characteristic for fluidic oscillators 7,32 . For the remaining sections that follow, all computations are conducted for Re = 350.

Phase reduction analysis.
In this part, we show how dimensionality reduction is achieved by employing phase reduction analysis, which reduces the theoretically infinite-dimensional multiphysics system into a single scalar phase equation. We first consider two configurations of the elastic structure placement. The first configuration has the elastic structures attached around the feedback channels and the second one has the structures attached around the mixing chamber, as shown in Fig. 4a, b respectively. For each condition, we use what is known as the direct method 24,33 to evaluate the phase sensitivity function Z of the system by applying a local- represents a gaussian impulse function in spatial coordinates with δ(t − t 0 ) being its temporal domain counterpart, and e y is a unit vector in the y-axis direction. The location x 0 in which these perturbations are applied, is shown in both Fig. 4a, b, and is imposed at multiple phase values determined by t 0 . Refer to the "Methods" section on how the phase sensitivity function is evaluated after the impulse perturbations are applied. The phase sensitivity functions have larger peak-to-peak values for the first configuration compared to the second one, as shown in Fig. 5a, b, respectively. This implies that the system is more sensitive to external perturbation applied around the feedback channels rather than those applied around the mixing chamber. Hence, we  www.nature.com/scientificreports/ can expect a wider synchronization region with respect to the amplitude and frequency of the external periodic forcing for the first configuration. In a more general context, these findings also show how placement of the external forcing around the feedback region of a system will cause a larger change in the system behavior. While the feedback region is easily identified for the case of the feedback-type fluidic oscillator, it is not always easy to identify the corresponding region for other kinds of oscillatory flows, such as for the case of vortex shedding 1,24 or feedback-free fluidic oscillator 34 . It is also found that there is a symmetrical property between the phase sensitivity functions of the top and bottom sides where Z top θ, − e y = Z bot θ + π, e y . This is expected given the spatially symmetric property of the oscillatory flow inside the fluidic oscillator and the phase difference introduced by the alternating sweeping motion. This finding is consistent with other studies of oscillatory flows with symmetrical properties 35 . This symmetrical property also implies that the largest synchronization region with respect to the external forcing frequency and amplitude is achieved when anti-symmetric forcing is applied to the system with a 1:1 synchronization. In other words, when the ratio of the external forcing frequency ω f and the original oscillation frequency ω n is approximately equal to 1. On the other hand, for symmetric forcing, it can be inferred that a larger synchronization region is likely to be achieved for with 1:2 synchronization ( ω f /ω n ≈ 2) 21,23 .

Synchronization analysis.
In this section, we validate the phase equation obtained from phase reduction analysis, and analyze the effects of periodic external perturbations on the oscillatory characteristics. From here onwards, we consider only the first configuration as shown in Fig. 4a. We analyze the phase dynamics of the system given that a periodic perturbation is applied in the form of a simple sinusoidal function. Utilizing the  www.nature.com/scientificreports/ symmetrical property, we apply anti-symmetric forcing to the outer boundaries of the elastic structures in the following form Using phase reduction theory, we can determine the synchronization region based on the external forcing frequency and amplitude, also known as Arnold Tongue 36 . From the results obtained with direct numerical simulations (DNS), synchronization is determined where the phase difference between the original oscillation and the external forcing converges to a constant value. To validate the synchronization region, multiple simulations were conducted using parameters 0.8 ≤ ω f /ω n ≤ 1.2 with step size of 0.02 , and 0.001 ≤ ε ≤ 0.005 with step size of 0.001 . The estimated synchronization region evaluated using phase reduction theory is in good agreement with results from DNS, as shown in Fig. 5c.
We further analyze the change in the physical characteristics of the oscillatory flow due to the periodic forcing. This is done by comparing the ratio of the root-mean-square (RMS) values and the maximum values of both C P and C D after the periodic forcing is applied with their nominal values when no external forcing is present. For the synchronizing cases, the RMS and maximum values can be evaluated after one period of oscillation once synchronization is observed. For the non-synchronizing cases, the phase difference will increase indefinitely 33,37 , so the said values are evaluated after at least one cycle of phase slip 33 is observed. The results are displayed in Fig. 6, with C P and C D denoting the ratio of RMS values and C max P and C max D denoting the ratio of maximum values. The parameters used are the same as the ones used for validating the synchronization region, and a smooth colormap is obtained with interpolation in-between the data points. For the synchronizing cases where ω f /ω n > 1 , it can be inferred that the overall internal fluid velocity also increases given that the oscillatory flow frequency increases. However, Fig. 6a, b, show that both the RMS and maximum values of C P decrease with increasing frequency ratio, which is somewhat contradictory with the conventional case where a larger outlet to inlet pressure difference is required to increase the flow velocity. We propose that this observation is obtained because the vibrating elastic , for the bottom side. www.nature.com/scientificreports/ structure imposes additional momentum into the fluid flow, such that a faster sweeping motion is maintained even without an increase in the inlet pressure. Although an increase in energy efficiency is seen in the fluid flow perspective, an additional energy source is required to actuate the elastic structure vibration. Even so, our proposed method avoids the need to increase the inlet pressure to increase the oscillation frequency. It is also worth noting that for non-synchronizing cases, the maximum value of C P increases. As mentioned, this can be undesirable since the material used to make the device needs to withstand a larger inlet pressure value 12 . This exemplifies the importance of estimating the synchronization region using phase reduction theory. We can also infer from the results that an increase of momentum in the lateral direction around the mixing chamber is observed for the synchronizing cases where ω f /ω n > 1 , as shown by the increase of both the RMS and maximum values of C D in Fig. 6c, d. The increase in oscillation frequency is a desirable effect for fluid mixing applications, since it is known that a higher oscillation frequency improves the mixing efficiency 7,38 . Since in the current study we consider an incompressible fluid such that its density does not change, the increase in lateral momentum indicates that there is an increase in the lateral fluid velocity, which promotes convective mixing 39 . We also note that the increase in the RMS value of C D is more apparent in the synchronizing cases compared to the non-synchronizing cases for the same value of periodic forcing amplitude where ω f /ω n > 1 . We can infer that a larger increase in the internal fluid velocity is achieved with a smaller forcing amplitude when synchronization is achieved. This further shows the advantage of estimating the synchronization region in advance and reduces the number of parametric studies required to find the optimal forcing amplitude and frequency to reach the desired system characteristics.
We have shown that by using phase reduction theory, the phase response of the given fluidic oscillator design is estimated well, and that synchronization with external forcing imposed by elastic structure vibration changes the internal fluid flow properties. This can be utilized in practical applications such as fluid mixing. Studies on mixing of chemical substances have considered the use of fluidic oscillators with many shapes and forms 7,11 . However, the resulting behavior is largely determined by the geometry of the device. The method proposed in this study introduces added controllability to the mixing behavior while retaining the overall geometry of the device. Moreover, the dimensionality reduction using phase reduction theory presented in this study might prove useful for other studies using fluidic oscillators, such as synchronization analysis of cascaded fluidic oscillators and controllability analysis by means of a different perturbation method 9,15,16 .
While several recent works have involved rigorous mathematical work to obtain analytical solutions for the phase equation of oscillatory flows 40,41 , we acknowledge that phase reduction alone does not allow estimation of the changes in the fluid flow properties due to synchronization. For instance, direct numerical simulations still need to be done for multiple values of periodic forcing frequencies and amplitudes before the changes in fluid flow properties can be analyzed, such as values of C P and C D as shown in Fig. 6. Phase-amplitude reduction 42,43 , which is a well-known extension of the phase reduction method, is a promising solution to this problem since it allows prediction of the future state variables such as the fluid flow properties. Although a sophisticated method of phase-amplitude reduction has been proposed recently 44 , direct implementation and confirmation of its validity for spatio-temporal dynamics in multiphysics system, including fluid-structure coupled dynamics, still needs to be confirmed.
Optimal control. The phase reduction analysis allows us to predict the synchronization condition. However, there is no constraint regarding the time required to reach the synchronization. Depending on the amplitude and shape of the periodic forcing signal, it can take multiple periods of oscillation before synchronization is achieved. This is why an optimal control framework 37,45,46 is desired with the goal to design a periodic forcing signal (i.e. the control signal) that achieves synchronization within one period of oscillation for a given averaged power.
The optimal control signal can be designed given that the phase sensitivity function Z , the phase difference during synchronization �ω , and its average power P are known, such that where η ε (ψ) is the optimal periodic forcing signal with amplitude ε , Z ′ (ψ) = dZ(ψ)/dψ is the gradient of the phase sensitivity function along the periodic phase value ψ = [0, 2π) , which can be evaluated numerically if Z is known, and is the Lagrange multiplier introduced for the optimization problem, and �ω = ω n − ω f . Explanation for Eq. (3) is described in the "Methods" section. The evaluated Z ′ (ψ) in this study is shown in Fig. 7a.
To demonstrate this idea, we choose a test case where a pure sinusoidal periodic perturbation with amplitude ε = 0.004 and frequency is applied to the system. When synchronization is achieved, the frequency difference is found at �ω = −0.06 and the calculated average power of the periodic forcing signal is P = 8 × 10 −6 . Given these values, we can construct the optimal control signal based on Eq. (3), as shown in Fig. 7b. The optimal control signal has a complex shape, since it contains many harmonic components. This is a consequence of the framework where η ε (ψ) is essentially a scaled version of Z ′ (ψ) , which contains many harmonic components in this case. However, judging from the amplitudes, each harmonic components exhibit only a small amount of power. Hence, we can adopt a reduced version of η ε (ψ) by decomposing it into its Fourier series components such that www.nature.com/scientificreports/ where a k and b k are the amplitude and phase shift of each k-th harmonic component respectively, and N represents the minimum number of harmonic components required to approximately represent the signal in Fourier series relative to its average power. For instance, we consider a constraint that 80% of the average power is required to reach optimal synchronization time. We can construct a control signal η ε (ψ) by choosing the M components of the Fourier series representation of η ε (ψ) that satisfies The resulting control signal after this reduction is shown in Fig. 7c, where the periodic forcing signal become less complex in shape, indicating that a practical implementation is more plausible.
We conduct direct numerical simulation to compare the time required for synchronization for two different cases, where the periodic forcing applied to the system is either a purely sinusoidal signal or the reduced version of the optimal control signal. Both periodic forcing satisfies �ω = −0.06 and P = 8 × 10 −6 as described previously in the design process of the optimal control signal. Figure 8a shows that it takes around 5 periods of oscillation to reach synchronization for a purely sinusoidal forcing, while Fig. 8b shows that the reduced version of the optimal control signal only requires one period of oscillation. We also note that there is no significant difference with the amplitude response of the driving force coefficient C D for both cases.
In general, when a dynamical system with a large degree-of-freedom is reduced to a phase equation, the resulting phase sensitivity function will have a complex waveform 47 due to the existence of higher-order harmonic terms. This also applies to the fluid-structure coupled system in the present study. Here, we proposed a method to optimize the transition time to the synchronized state by designing a smoother waveform of the external forcing and demonstrated its effectiveness. The proposed method can reach a synchronized state with optimal transition time, even for actual experimental systems where a complex periodic waveform cannot be applied due to physical constraints.

Conclusions
In this work, we have introduced a design concept of a feedback-type fluidic oscillator surrounded by elastic structures. Our approach maintains the intrinsic oscillatory flow characteristics inside the fluidic oscillator but introduces controllability through synchronization with weak external forcing from the moving component. The resulting changes in physical characteristics due to synchronization can be advantageous for practical applications such as fluid mixing, without any need to change the overall geometry or to increase the inlet pressure. We demonstrate how the synchronization condition can be estimated well using phase reduction theory, even for a complex fluid-structure coupled system, allowing for the desired system characteristics to be discovered more efficiently. Furthermore, taking practical applications into account, we demonstrate how synchronization can be achieved within one period of oscillation by implementation of the optimal control framework. The geometry of the computational model is shown in Fig. 2 where its dimensions are normalized to the width of the feedback channel d . The power nozzle width is d i = 0.7d . The elastic structure on each side has a thickness of 0.2d . A parabolic velocity profile U in = 6U(y + 1) H − (y + 1) /H 2 is set for the fluid inlet boundary condition, where U is the mean velocity and H = 2d is the inlet channel width. The fluid outlet boundary condition is set to a static pressure p out = 0 . The left and right boundaries of the elastic structures are fixed, that is ∂w s /∂t = 0 . The computational domain is filled with quadrilateral meshes that are symmetrical with respect to the line y = 0.

Methods
Phase reduction and synchronization analysis. We briefly describe the formulation for phase reduction analysis as follows. Consider the system state of the FSI dynamics X(x, t) = [u f (x, t), w s (x, t)] T and the state dynamics when an external forcing is applied as where x = x, y is the coordinates in which measurements are taken that is theoretically infinite, F{·} is a functional representing the state dynamics, ε is the perturbation amplitude and η(x, t) is the perturbation function. The system exhibits a stable limit-cycle oscillation when no perturbation is applied (ε = 0) , such that X(x, t) = X(x, t + 2π/ω n ) . By introducing a phase functional {X(x, t)} , the system state can be mapped into scalar phase values such that θ(t) = {X(x, t)} , where θ ∈ [0, 2π] . The phase functional is defined in such a way that θ increases with a constant frequency ω n , that is where Z(θ; x) = δ /δX is the phase sensitivity function, which can be approximated along the limit-cycle orbit. Should an external forcing with amplitude 0 < ε ≪ 1 be applied, the phase dynamics now becomes and ∇ · u f = 0, www.nature.com/scientificreports/ Given this formulation, we can utilize impulse perturbation method, also known as the direct method 24,33 , in order to evaluate Z(θ; x) . We consider the case in which a weak and localized impulse forcing is applied, that is ε ≪ 1 and η(x, t)dx = δ(t − t 0 ) , where δ is a gaussian impulse function with small enough standard deviation. In this study, we achieve this by setting η(x, t) = δ(x − x 0 )δ(t − t 0 ) − e y . The gaussian impulse function δ(x − x 0 ) has a standard deviation of 0.1d , while δ(t − t 0 ) has a standard deviation of 0.02T , where T is the oscillation period for C D . By applying the impulse perturbation, a phase shift asymptotic to the original limit cycle orbit g(θ; x, ε) will be introduced. We can then evaluate the phase sensitivity function by the following approximation.
Moving on, we will consider the case where the perturbation is localized in the spatial coordinate, such that x can be dropped from the notation. We applied the impulse perturbation at 20 different phase values within a single oscillation (determined by t 0 ) in order to obtain the phase sensitivity function Z , as shown in Fig. 5a, b. The synchronization condition is achieved if the phase difference between the oscillatory system and the periodic forcing φ(t) = θ(t) − ψ(t) converges to a constant value (phase-locked), where ψ(t) = ω f t . In other words, dφ/dt = �ω + εŴ(φ) converges to 0, where �ω = ω n − ω f and Ŵ(φ) is the phase coupling function defined as which is a 2π periodic function that is evaluated using averaging approximation due to the fact that φ(t) is changing slowly over time compared to ψ(t) . With these formulations, it is possible to approximate the synchronization condition based on the forcing amplitude and frequency, such that Optimal control framework. Recall that synchronization is achieved when the oscillatory system and the periodic forcing are phase-locked. By this definition, a stable fixed point φ * s must exist, such that dφ/dt = �ω + Ŵ ε (φ) converges to 0 at φ = φ * s , where Ŵ ε (φ) = εŴ(φ) for simplicity. We can infer that in order to achieve the fastest convergence, we need to have the steepest gradient possible for Ŵ ε (φ) around φ = φ * s . In other words, the goal of optimal control is to This goal must be combined with the following phase-locking constraint and the averaged power constraint where P is the averaged power exhibited by the periodic forcing signal η ε (ψ) with amplitude ε , and �·� denotes the averaging operator. Through these goals and constraints, we can construct the following cost function of the optimization problem Furthermore, by using calculus of variations, the optimal periodic forcing signal can be proven to be equal to Eq. (3). Equation (3) also shows that the optimal periodic forcing signal can be designed given that the phase sensitivity function Z , the phase difference during synchronization �ω , and its averaged power P are known. In practice, we can initially figure out Z using the phase reduction analysis and use a purely sinusoidal periodic forcing signal to figure out P which is required to reach synchronization with a particular �ω based on the synchronization condition shown in Eq. (14). These values can then be substituted back to Eq. (3) to evaluate the optimal periodic forcing signal η ε (ψ).

Data availability
The experimental dataset that supports the findings of this study is available from the corresponding author upon reasonable request.