Modulation of viscoelastic fluid response to external body force

Transient flow responses of viscoelastic fluids to different external body forces are studied. As a non-Newtonian fluid, the viscoelastic fluid exhibits significant elastic response which does not raise in Newtonian fluid. Here, we investigate the transient response of a viscoelastic Poiseuille flow in a two-dimensional channel driven by external body forces in different forms. The velocity response is derived using the Oldroyd-B constitutive model in OpenFOAM. Responses in various forms like damped harmonic oscillation and periodic oscillation are induced and modulated depending on the fluid intrinsic properties like the viscosity and the elasticity. The external body forces like constant force, step force and square wave force are applied at the inlet of the channel. Through both time domain and frequency domain analysis on the fluid velocity response, it is revealed that the oscillation damping originates from the fluid viscosity while the oscillation frequency is dependent on the fluid elasticity. The velocity response of the applied square waves with different periods shows more flexible modulation signal types than constant force and step force. An innovative way is also developed to characterize the relaxation time of the viscoelastic fluid by modulating the frequency of the square wave force.

Fluid has been widely used as the medium for energy and force transmission in control systems named as fluidic devices. Such devices has been developed into functional dynamic elements like amplifiers and triodes, and exhibit logical operations similar as electronic components [1][2][3][4][5] . Much attention has been paid to the fluid logic component due to its advantages of no extra mechanical structure, low cost and high reliability. Its resistance to extreme conditions of shock, vibration or radiation is also desired for many applications.
Traditionally, fluidic logical functions are implemented in macro-scale devices at high Reynolds (Re) number flows by taking the advantage of the nonlinear property of Newtonian fluid 6 . While in microfluidic component, the logical function usually relies on the formation of droplets 7 , bubbles 8 or circuits with valve switches 9 . For real application, these methods usually require complex control systems. On the other hand, realizing nonlinear behavior of Newtonian fluid for microfluidic components remains large challenging as extremely high pressure is required to achieve high Re number condition 10 . However, other studies [11][12][13][14] show that viscoelastic fluid can induce nonlinear characteristic from its elastic property even at extremely low Re conditions [15][16][17] . Therefore, by using a viscoelastic fluid, a nonlinear microfluidic circuit can be obtained.
Based on the viscoelastic fluid, many microfluidic logical components are developed for microfluidics applications [18][19][20] . Microfluidic diodes, rectifiers and amplifier are three most common microfluidic circuit components functioning with viscoelastic fluid [20][21][22] . Grosiman et al. designed a bistable flip-flop memory and a flux stabilizer with viscoelastic fluid in a continuously curved sudden contraction structure 21 . Also, they proposed a microfluidic rectifier device 22 using a nozzle and diffuser structure by taking the viscoelastic flow direction as the logic operation signal. It had been successfully used to characterize the viscoelastic rectifier [23][24][25] . The above studies mainly took the viscoelastic flow direction as the logic operation signal without investigating the response of fluid itself. Actually, the transient response of the viscoelastic fluid to external excitations will provide a new way to implement reliable microfluidic circuit components with dynamic modulation functions. The response of viscoelastic fluid has been previously studied for a start-up flow under an impulse or constant pressure gradient at rest [26][27][28][29] . The start-up flow for non-Newtonian fluid are particularly for the verification of the numerical methods, which are used for calculating transient flow responses [30][31][32] .

Formulation and Numerical Method
Physical model. This paper studies the velocity response of viscoelastic fluid to an applied external body force in a two-dimensional fully-developed Poiseuille flow model as shown in Fig. 1. The widths of both the inlet and the outlet 2H are set at a dimensionless number of 2; the flow distance L is set at 8. The two walls are in the non-slip boundary condition. An external body force is exerted in the x-direction.
Governing equations for viscoelastic fluid flow. Firstly, we need to choose suitable constitutive models to represent the flow properties of viscoelastic fluids. There are various constitutive models being developed, among which two models are often used 34 36 , etc. In this study, we adopt the Oldroyd-B model, which is one of the Oldroyd 8-Constant models due to its simplicity. Based on the Oldroyd-B model, we build a new solver for viscoelastic fluid flow in OpenFOAM.
The dimensionless governing equations for incompressible viscoelastic fluid flow are expressed as:: where u is the velocity vector normalized by the maximum velocity U max ; p is the pressure in the fluid and normalized by ρU max 2 , where ρ is density of the fluid; T s is the viscous stress, and T p is the elastic stress. T s is described as where C is the conformation tensor; I (Kronecker symbol) is the unit tensor; η P is the dynamic viscosity of the solute, and = λ Wi u L is the Weissenberg number which describes the strength of the elasticity of the fluid, λ is the relaxation time of the viscoelastic fluid. The elastic stress term is a symmetric positive definite second-order tensor.
The conformation tensor transport equation is:  www.nature.com/scientificreports www.nature.com/scientificreports/ Eq. (4) the hyperbolic characteristic of Eq. (4), an iterative algorithm, and the existence of equilibrium. The rapid growth rate of conformation tensor is inexact to polynomial fitting, resulting in divergence of calculation 37,38 . Instead of solving the elastic stress tensor or the conformation tensor, we solve the log-conformation tensor, where the extensional components of deformation field act additively rather than multiplicatively. After solving the log-conformation tensor, the conformation and elastic stress tensors can be recovered.
In this paper, we reconstructed the conformation tensor by the alternative exponential method, an algorithm that turns the nonlinear terms in the equations into linear terms. Specifically: T The velocity gradient can be expressed as Ω ∇ = + + − u B NC 1 , which decomposes velocity gradient into extensional component B and rotational components Ω and N. The C, B, Ω, N, and D are described as follows: The conformation tensor transport equation can be transformed as follows: Here, we introduce a parameter Ψ associated with R and Λ, where R is the feature vector of Ψ, Λ is diagonal matrix constituted by the eigenvalues of Ψ.
T Substituting Eqs (9) and (10) into Eq. (7), we get: The dominant nonlinear term of original conformation tensor in Eq. (4) turns into a linear term through strain rate tensor reconstruction in Eq. (11), which can significantly improve the computational stability. LCR-based solver performance has been verified in our previous study 39 . Numerical procedure. The numerical process is shown in Fig. 2. In the numerical calculation, an implicit finite-volume approach is used; the first-order Euler scheme is applied for time marching of all unsteady equations 40 with the dimensionless time step δt being 10 −3 . In spatial discretization, the QUICK scheme 41 is adopted for the convective term in Eq. (2), and the MINMOD scheme 42 are used for the logarithmic conformation tensor transport. The relationship between different external forces and the velocity response signal can be obtained from one-dimensional flow by comparative analysis.

Results and Discussion
In this paper, we investigate the velocity response of viscoelastic fluid in the two-dimensional fully-developed Poiseuille flow under external forces in different forms, i.e., a constant force, a rectangular function force, and a square wave force.

Velocity response of viscoelastic fluid flow under a constant force.
In this section, we first compare the transient velocity responses between the Newtonian fluid and the viscoelastic fluid flows when a constant force is applied. We set the applied constant force F 0 = 2 from time t = 0 and β is 0.01. The elasticity number, which is defined as the ratio of elastic stress to the inertial stress, or E = Wi/Re, is set as E N = 0 for the Newtonian fluid and E V = 1 for the viscoelastic fluid.
The centerline velocities u (dimensionless) of two fluid flows at the outlet in a dimensionless time t are shown in Fig. 3. For the Newtonian fluid flow, the velocity u N first increases logarithmically from time 0 to 1 and then saturates at u N = 1 due to the balance between the applied force and the viscous dissipation in the flow. In contrast, www.nature.com/scientificreports www.nature.com/scientificreports/ the velocity of viscoelastic fluid u V fluctuates as an under-damped vibration function with time: u V first increases to a peak value of 1.9 at t = 1.5, drops to a dip value of 0.7 at t = 4.1, and then continues to oscillate with a damping effect until it finally stabilizes at u V = 1. The oscillation originates from the non-zero elasticity of viscoelastic fluid responding to the constant F 0 . u V increases at time 0 due to the applied F 0 . At the same time, the molecular microstructures in viscoelastic fluid are deformed due to its elasticity and elastic potential energy is stored in the molecular microstructures of viscoelastic fluid. The potential energy can be converted to kinetic energy which makes u V continue to increase and exceed the stable velocity of the Newtonian fluid until time t = 1.5. And then, u V starts to drop due to the viscous effect. Meanwhile, the kinetic energy of fluid flow converts to the potential energy, which makes u V continue to drop below the stable velocity of the Newtonian fluid flow until t = 4.1. At this moment, F 0 overwhelms the viscous dissipation and u V starts to increase again. The process repeats which leads to an oscillation velocity response of viscoelastic fluid flow to the constant F 0 . During the oscillation, the total energy decreases due to the viscosity, and therefore the oscillation is damping. The damped oscillation of velocity response in viscoelastic fluid flow reaches a stable state when the viscous resistance balances with the external F 0 and the molecular microstructure in viscoelastic fluid remains in a static state. Then, viscoelastic fluid flows at a constant velocity u V = 1 after t = 15. The stable velocities are the same for the Newtonian fluid and viscoelastic fluid flows with the same viscosity under the same constant force. Further investigation on the damped oscillation in the frequency domain reveals a peak resonance at f = 0.27 based on fast Fourier transform, as shown in the inset of Fig. 3.
When the constant forces F 0 with different magnitudes are applied to viscoelastic fluid flows, the velocity oscillates with an amplitude proportional to F 0 while the oscillation peaks and dips occur at the same time value for all F 0 , as indicated by open circles in Fig. 4(a). Here, the oscillation is assumed to be a damped harmonic wave function and is expressed as t c 0 By fitting parameters in Eq. (12) to the numerically calculated values for different F 0 , the oscillation curves of u are shown as solid lines in Fig. 4(a). The fitted curves closely match with numerical results, confirming that the  www.nature.com/scientificreports www.nature.com/scientificreports/ velocity oscillates in a damped manner. The parameters are derived for different F 0 as shown in Table 1, where A and T are the amplitude and period of harmonic oscillation, respectively; α is the damping coefficient; u 0 and t c are the offset values for the amplitude and time, respectively. α and T remain the same when F 0 changes because they are intrinsic properties of viscoelastic fluid and only affected by the viscosity and elasticity of viscoelastic fluid. The amplitude of sinusoid functions A, on the other hand, increases proportionally with F 0 , which is plotted in Fig. 4(b). The time offset t c , which relates to the oscillation phase of velocity, is the same for different F 0 because the starting time of force is the same.
As analyzed above, the damped oscillation depends strongly on the elasticity of viscoelastic fluid. In a fluid with larger elasticity, more potential energy can be stored in the fluid before it is saturated. The potential energy is then converted to the kinetic energy and drives the velocity to increase for much time. For the same reason, with larger elasticity, more time is also needed for the fluid to convert the kinetic energy to the potential energy during the velocity drop period. As a result, the oscillation period should be longer for the fluid with a larger elasticity number. This is numerically confirmed by comparing the transient response of centerline velocities at different elasticity numbers as shown in Fig. 5(a). The colour bar represents the velocity amplitude. It is found that the fluid with larger elasticity oscillates slower than the fluid with less elasticity. Also, the oscillation amplitude is larger for the fluid with larger elasticity because more energy can be stored in this fluid and be converted to the kinetic energy.
The velocity response at different elasticity numbers is then fitted using Eq. (12) and the coefficients are shown in Table 2. The oscillation period T increases from 4.16 to 17.94 as elasticity number E increases from 1 to 20. This proportional relation is explained quantitatively previously. The amplitude A also increases significantly from 1.79 to 7.26, which is because of a viscoelastic fluid with larger elasticity gains and converts more potential energy to the kinetic energy. In addition, the effective damping factor T· α, or the damping in one oscillation period, decreases from 2.20 to 0.66, indicating that the fluid with larger elasticity will experience a relatively lower damping effect under a constant force excitation.
The oscillation is further investigated in the frequency domain through a fast Fourier transform of temporal response. Different resonant frequencies of f = 0.242, 0.110, 0.078, 0.064 and 0.055 are obtained for E = 1, 5, 10, 15 and 20, respectively, as shown in Fig. 5(b). We refer to this resonant frequency as the dominant frequency of viscoelastic fluid. The dominant frequency drops as the elasticity number increases, as seen in Fig. 5(c), which results from the increase of oscillation period. Meanwhile, the peak amplitude of velocity increases with the elasticity number, as shown in Fig. 5(d). This trend is caused by the more kinetic energy converted to the potential energy during the molecular deformation process and more potential energy converted to the kinetic energy during the molecular release process.  www.nature.com/scientificreports www.nature.com/scientificreports/ Besides elasticity, the viscosity is another critical property of viscoelastic fluids which describes its resistance to gradual deformation under shear stress or tensile stress. This resistance induces friction, which slows down the flow velocity. The centerline velocities are compared for viscoelastic fluids with a fixed elasticity number E = 1 and different β in Fig. 6(a). As expected, the oscillation amplitude of velocity drops faster in viscoelastic fluid with larger viscosity than that in one with lower viscosity. It can be also observed that the peaks and dips of oscillation occur at the same time interval regardless of viscosity. Therefore, the oscillation period does not change with the viscosity, but the damping behaviour does. The fitting coefficients obtained using Eq. (12) for different β are listed in Table 3. The oscillation period T and the amplitude A remain almost the same when β increases from 0.01 to 0.2. In contrast, the damping factor α increases significantly from 0.06 to 0.29 for the same values of β, suggesting that the oscillation damping depends remarkably on the viscosity. The dominant frequency is the same for the viscoelastic fluid with different β as shown in Fig. 6(b) because their oscillation period is not affected, however, the amplitude at the dominant frequency is the highest for the fluid with the lowest β, which is due to the lowest energy loss in the oscillation.   www.nature.com/scientificreports www.nature.com/scientificreports/ Velocity response of viscoelastic fluid flow under a rectangular function force. In this section, we study the transient velocity response of viscoelastic fluid flow with elasticity E = 1 to an applied rectangular function force. Here, we define the force as: r r r 0 where T r is the force applied time. Figure 7(a) shows the transient velocity response to the rectangular force with different F 0 and fixed T r = 1. The velocity increases from time t = 0 to t = T r , and during the period, the condition    www.nature.com/scientificreports www.nature.com/scientificreports/ is the same with one when a constant force is applied. As discussed previously, under a constant force, the velocity oscillation for the fluid with E = 1 has a period T = 13.7 and is in its ramp-up stage of the first cycle when t = T r . However, for the rectangular force, the force is unloaded at this point. One interesting phenomenon observed is that the velocity remains constant for another 1 s after the force is unloaded until t = 2. After that, the velocity amplitude starts to drop sharply from t = 2 to t = 4 and then drop very slowly from t = 4 to t = 6. Then it becomes damped oscillation and finally ends with a velocity of 0. The length of this "platform effect", the short-term constant velocity behaviour after the force is unloaded, is independent of the rectangular force amplitude but the magnitude of velocity increases proportionally with the amplitude. A fast Fourier transform is performed to further investigate the velocity response, as shown in Fig. 7(b). Responses to rectangular forces of all amplitudes exhibit a resonance frequency f = 0.091 Hz, which is blue shifted compared to the resonance frequency f = 0.078 Hz under the constant force excitation condition as indicated in Fig. 6.
By using viscoelastic fluid with a larger elasticity (E = 10) and varying T r , it continues to study the platform effect. It is found that the platform effect only appears when the force is unloaded during the velocity ramping stage of the first oscillation period, as seen in Fig. 8(a). At this stage, the molecular microstructure of viscoelastic fluid is stretched under large shear stress if an external force is applied and potential energy is stored in viscoelastic fluid. Due to the memory characteristic of viscoelastic fluid, the potential energy releases, and the velocity is maintained at the same value until the conversion from the potential energy can no longer support the flow at this velocity, which will, therefore, start to drop afterwards. The energy conversion time is related to the elasticity  www.nature.com/scientificreports www.nature.com/scientificreports/ of viscoelastic fluid, which ends at the ¼ period, or the peak value time of the first oscillation. Therefore, it can be concluded that a viscoelastic fluid with larger elasticity will display a longer platform effect due to the longer period. In the fast Fourier transform of the response, resonance occurs at the same frequency value of f = 0.091 Hz for different T r , which indicates that the domain frequency does not change with T r . In contrast, the power spectra density at the resonant frequency changes significantly when the unloading time changes. Therefore, the rectangular force does not affect the domain frequency but only determines the resonant amplitude.
Velocity response of viscoelastic fluid flow under a square wave force. In this section, a square wave force F s is applied to the viscoelastic fluid, which is defined as where T s is the period of the square wave. We investigate the velocity response to the square wave force on a viscoelastic fluid with elasticity E = 1. Here, square wave forces with period T s = 1, 4, 8 and 16 are applied to the viscoelastic fluid, and the velocity response is shown in Fig. 9. For T s = 1, 8 and 16, multiple resonance peaks are observed. These peaks deviate greatly from the intrinsic period of the velocity oscillation T i = 4.16, which is found earlier for a viscoelastic fluid under a constant force (as shown in Table 2). At T s = 4, which is close to T i , each period has a single peak, which is consistent with single resonance. We further study the response of viscoelastic fluid flow to a square wave force in the frequency domain. A Fourier transform is applied to each velocity response for square waves with periods of 1, 4, 8 and 16 s. The applied square wave excitation is an odd function and its corresponding Fourier expansion can be written as 39 where f = 1/T s . Therefore, resonant peaks are induced at n·f, where n is equal to 1, 3, 5 and so on. For example, resonant peaks appear at 1, 3, and 5 Hz when T s = 1 with a corresponding f of 1 Hz as shown in Fig. 10(a). Similarly, for T s = 8 and the corresponding f = 0.125 Hz, resonant peaks appear at f (0.12 Hz), 3 f (0.38 Hz) and 5 f (0.63 Hz), as shown in Fig. 10(c). By the same mechanism, at T s = 16 s or f ≈0.06 Hz, resonant peaks are induced at f (0.06 Hz), 3 f (0.19 Hz) and 5 f (0.31 Hz) and so on. However, when T s is close to the natural period of viscoelastic fluid, or 4 s, other order frequencies are significantly suppressed, and only the first resonant mode at f = 0.24 Hz is strongly induced, as shown in Fig. 10(b). This is because the stimulus frequency matches the natural frequency and most of the energy is coupled at this frequency. As a result, the resonant amplitude is also much larger than that in other cases. Because the domain frequency of viscoelastic fluid flow can be tuned just by changing the corresponding elasticity, such fluids can be readily used for the narrowband filtering and amplification of signals. This is especially useful in biochemical applications due to fluid compatibility.

Conclusions
We investigate the transient velocity response of viscoelastic fluid flow to externally applied forces based on numerical simulations. The velocity responses to different forces are interpreted which highly depend on the intrinsic properties of viscoelastic fluid, particularly its elasticity and viscosity. Some interesting results are obtained as follows: (1). A damped oscillation velocity response to a constant force is obtained. The oscillation amplitude, period and damping coefficient can be modulated by changing the viscosity coefficient β and the elasticity E of the viscoelastic fluid. It demonstrates that the period of the oscillation signal relies on the elasticity number of the fluid while the damping coefficient is dependent on viscosity. (2). We investigate the velocity response to a rectangular force and a platform effect in a viscoelastic fluid was observed, by which the velocity is maintained at a constant value for certain time after the force is unloaded. (3). A periodic square wave force is also applied and harmonic resonance is occurred when the applied force period matches the intrinsic period of the viscoelastic fluid. The resonance is due to the coupling between the force and the fluid and the resonant amplitude of the velocity is maximized.
By analyzing the response features of viscoelastic fluid excited by external forces, the transient velocity can be utilized as modulated signal when processing different force functions in microfluidic circuit. In general, the nonlinear response to the external force allows the viscoelastic fluid to be widely applied in the waveform modulation and filter applications, such as digital signal processing, signal converter, controller, signal encryption, password converter, etc.