Numerical Investigation of T-Shaped Microfluidic Oscillator with Viscoelastic Fluid

Oscillatory flow has many applications in micro-scaled devices. The methods of realizing microfluidic oscillators reported so far are typically based on the impinging-jet and Coanda effect, which usually require the flow Reynolds number to be at least at the order of unity. Another approach is to introduce elastomeric membrane into the microfluidic units; however, the manufacturing process is relatively complex, and the membrane will become soft after long-time operation, which leads to deviation from the design condition. From the perspective of the core requirement of a microfluidic circuit, i.e., nonlinearity, the oscillatory microfluidic flow can be realized via the nonlinear characteristics of viscoelastic fluid flow. In this paper, the flow characteristics of viscoelastic fluid (Boger-type) in a T-shaped channel and its modified structures are studied by two-dimensional direct numerical simulation (DNS). The main results obtained from the DNS study are as follows: (1) Both Weissenberg (Wi) number and viscosity ratio need to be within a certain range to achieve a periodic oscillating performance; (2) With the presence of the dynamic evolution of the pair of vortices in the upstream near the intersection, the oscillation intensity increases as the elasticity-dominated area in the junction enlarges; (3) Considering the simplicity of the T-type channel as a potential oscillator, the improved structure should have a groove carved toward the entrance near the upper wall. The maximum oscillation intensity measured by the standard deviation of flow rate at outlet is increased by 129% compared with that of the original standard T-shaped channel under the same condition. To sum up, with Wi number and viscosity ratio within a certain range, the regular periodic oscillation characteristics of Oldroyd-B type viscoelastic fluid flow in standard T-shaped and its modified channels can be obtained. This structure can serve as a passive microfluidic oscillator with great potential value at an extremely low Reynolds number, which has the advantages of simplicity, no moving parts and fan-out of two.


Introduction
A fluid oscillator is a device that utilizes the instability of fluid flow to produce a continuous pulsating jet. It has a wide variety of applications, such as jet mixing enhancement [1][2][3], heat transfer enhancement [4][5][6], cavity resonance suppression [7,8], etc. In addition, as a flow control actuator, the fluid oscillator is utilized in aerospace industries because of its unique features and its simplicity with no moving parts, such as missile control [9] and jet thrust vector control [10]. A comprehensive overview of the fluid oscillator can be found in the review literature [9,11].
The concept of the microfluidic oscillator has emerged since the size of the fluid oscillator device was reduced to microscale. There are many applications of microfluidic oscillator in microscale fields, such as particle sorting [12]/focusing [13], mixing/heat transfer enhancement at microscale [14,15], flowmeter [16][17][18], chemical process enhancement [19,20], filtration improvement [21] and hemodynamics [22]. Microfluidic oscillators can also be integrated with other fluidic components for more complex flow control and operation.
At present, the principle of using continuous fluid to realize the fluidic oscillator can be roughly divided into the following categories (as exemplified in detail by several typical literature articles in Table 1). (1) The equivalence principle of electronic and fluidic circuit. Lammerink et al. [23] designed several equivalent basic fluidic elements and assembled them into a multivibrator. (2) Coanda effect. A rapid jet enters the diffuser at a high aspect ratio and takes the nearby fluid element away from the wall. Therefore, the uneven low-pressure area is formed on both sides, and then the flow is randomly attached to one of the side walls. Xie et al. [24] and Yang et al. [25] studied the Coanda effect with symmetric and asymmetric feedback channels, respectively. (3) Impinging-jet-based principle. The flow instability can be induced by the jet impinging on a certain barrier [16,18] or by the two jets colliding at a certain angle [26,27], so as to realize the oscillating characteristics. The common features of the above design principles are that Newtonian fluid is used as the flow working fluid, and the flow Reynolds number (Re) ranges from about 10 to several thousands. However, as Re further decreases to an order of even smaller than unity, the inertial nonlinearity effect is usually so weak as to restrict the fluid flow to be laminar.
The core of a fluidic circuit is the nonlinearity [28]. As one of the components of a fluidic circuit at microscale, the microfluidic oscillator needs other non-inertial nonlinearity sources, which can be introduced in many ways. The methods of multiphase flow, such as micro bubbles [29,30]/droplets [31] and electrochemistry [32] for microfluidic control are complex in channel fabrication and lack stability or repeatability. Recently, a hydroelastic microfluidic oscillator [33,34], which uses an elastic diaphragm (silicone rubber) as the intermediate layer of the channel, was realized within a medium range of Re (~10-100). Similarly, a microfluidic oscillator was achieved from a channel with two valves made from an elastomeric membrane [35,36]. Nevertheless, the stiffness of the elastic diaphragm will decrease after a long period of operation, resulting in the overall reduction of the working pressure range of the hydroelastic microfluidic oscillator. In addition, the complexity of the channel increases the difficulty of its fabrication.
Despite the almost negligible inertial effects in a micro-scaled channel, the viscoelastic fluid can bring about complex flow phenomena, such as symmetry breaking, purely elastic instabilities [37,38] and even elastic turbulence [39,40], which arise from the combination of curved streamlines and large tensile stress. The anisotropic elastic stress arising from the natural nonlinear property of viscoelastic fluid flow itself is the source of strong nonlinearity, which meets the core requirements of fluidic circuit [28]. Taking the viscoelastic fluid as the working medium has the advantages of a single-phase medium and no moving parts compared with the abovementioned methods.
The application of oscillatory flow of viscoelastic fluid has been reported in several studies. Using the flow of viscoelastic fluid in the form of oscillation generated by external solenoid valve, Asghari et al. [41] achieved particle focusing and separation at nano-scale, which indicates that the same efficiency can be obtained by reducing the channel length from 4 cm to 4 mm. Sun et al. [42] studied the weakly shear-thinning fluid flowing in the microfluidic oscillator previously designed for Newtonian fluid. The results verified that the linear relationship between flow rate and frequency is still valid. All these studies involved viscoelastic fluid, but the typical properties of viscoelastic fluid have not been fully used to realize the oscillator.
The realization of a microfluidic oscillator via a simple structure is of great importance for eliminating the dependence on external instruments, integrating microfluidic components and reducing the failure rate of functional components. Therefore, the base geometry studied in this paper is a right-angled axisymmetric T-junction, which is a classical geomet-ric model used in hemodynamics [43], droplet formation [44], oil transportation [45], etc. The purpose of choosing symmetrical geometry is to expect the average flow rate of both outlets to be similar, so that both exits can serve as the excitation source of pulsation. In fact, the study of unsteady flow through a T-junction is fairly common. Soulages et al. [46] investigated the T-shaped stagnation flow geometries with and without a short length of recirculating cavity, and found the onset of the symmetry-breaking transition with the presence of the cavity. Varshney et al. [47] experimentally studied a millimeter-sized T-type channel with a long recirculating cavity. They reported that in the range of similar Re and Wi numbers where forward bifurcation usually occurs, forward Hopf bifurcation was observed instead, and the existence of the long recirculating cavity played a decisive role. Poole et al. numerically investigated the effects of aspect ratio and shear-thinning properties for generalized Newtonian fluid [48] and viscoelasticity [49] on the symmetrybreaking bifurcations in T-channel flows. A common feature of the above studies is that two streams collided head-on from two opposite entrances and then evacuated through 90 degrees. In terms of the two possibilities of flow in or out of each port in the T channel, two streams can collide in perpendicular [44] or one stream enters from the left and leaves from both the exit in line with the inlet and the exit perpendicular to the inflow [43,50,51], as summarized in Table 2. So, naturally, we were very curious about the last case, that is, the fluid flowing out from two outlets in a straight line. Although seldom, studies like Colin et al.'s [45] have also adopted this kind of arrangement; they have focused on the impact of wormlike micelles on the T-junction flow, which is a classical model in the oil industry. The elastic instability of the viscoelastic fluid flow in the T channel with two streams leaving in opposite directions is still worthy of further study.
In this paper, two-dimensional direct numerical simulation (DNS) of viscoelastic fluid flow in a T-shaped channel as prototype and several modified structures is carried out to investigate the flow characteristics at a very low Re. The effect of viscosity ratio of the viscoelastic fluid is also studied. The remainder of this paper is arranged as follows. Section 2 introduces the numerical simulation methodology, including the description of governing equations, the verification of grid independence, etc. In Section 3, the results and discussions are presented mainly from the pulsing performance, the effect of Wi and viscosity ratio, and the stretching analysis of macromolecular structure. Finally, Section 4 ends with the conclusions, and the prospect of the future improvement of a T-shaped microfluidic oscillator is put forward.         2D n/a Carreau-Yasuda 50~1000 n/a n/a [51] Two streams collide in perpendicular Exp 100 200, 600 and 1000 ppm CTAC with NaSal 1.54~42.22 n/a n/a [44] Two streams leave in opposite directions 3D 1000 modified Giesekus n/a n/a ~0 [45] Note: PEO: polyethylene oxide; CTAC: cetyltrimethylammonium chloride; NaSal: sodium salicylate; PAAm: polyacrylamide; MW: Molecular weight. a The abbreviations for Exp and Num denote the results of experiments and simulations, respectively. Simulation includes two types: two-dimensional (2D) and three-dimensional (3D). b The data is estimated from the raw data in the original paper. 2D n/a Carreau-Yasuda 50~1000 n/a n/a [51] Two streams collide in perpendicular Exp Micromachines 2021, 12, x FOR PEER REVIEW 8 of 26 2D n/a Carreau-Yasuda 50~1000 n/a n/a [51] Two streams collide in perpendicular Exp 100 200, 600 and 1000 ppm CTAC with NaSal 1.54~42.22 n/a n/a [44] Two streams leave in opposite directions 3D 1000 modified Giesekus n/a n/a ~0 [45] Note: PEO: polyethylene oxide; CTAC: cetyltrimethylammonium chloride; NaSal: sodium salicylate; PAAm: polyacrylamide; MW: Molecular weight. a The abbreviations for Exp and Num denote the results of experiments and simulations, respectively. Simulation includes two types: two-dimensional (2D) and three-dimensional (3D). b The data is estimated from the raw data in the original paper. 2D n/a Carreau-Yasuda 50~1000 n/a n/a [51] Two streams collide in perpendicular Exp 100 200, 600 and 1000 ppm CTAC with NaSal 1.54~42.22 n/a n/a [44] Two streams leave in opposite directions 3D 1000 modified Giesekus n/a n/a ~0 [45] Note: PEO: polyethylene oxide; CTAC: cetyltrimethylammonium chloride; NaSal: sodium salicylate; PAAm: polyacrylamide; MW: Molecular weight. a The abbreviations for Exp and Num denote the results of experiments and simulations, respectively. Simulation includes two types: two-dimensional (2D) and three-dimensional (3D). b The data is estimated from the raw data in the original paper. Figure 1 shows the schematic of the flow unit studied in this paper. The basic structure is a T-shaped channel with one inlet and two outlets. Four modified structures were formed by adding cavities. The geometry description and symbol marking of the channels are summarized in Table 3.

Physical Model
Micromachines 2021, 12, x FOR PEER REVIEW 9 of 26 Figure 1 shows the schematic of the flow unit studied in this paper. The basic structure is a T-shaped channel with one inlet and two outlets. Four modified structures were formed by adding cavities. The geometry description and symbol marking of the channels are summarized in Table 3.

Governing Equations and Numerical Method
Two-dimensional incompressible adiabatic flow was considered. The equations of conservations of mass and momentum were adopted. Generally, viscoelastic fluid has shear-thinning physical property, but to simplify the research complexity, Boger-type fluid with constant shear viscosity was selected in this paper, and the constitutive equation was Oldroyd-B [53]. This simplification is quite common for the convenience of research. [54,55]. For the sake of generality, each variable is nondimensionalized, and the dimensionless governing equations are the following.
where is the time, is the velocity, is the pressure, ∇ is the nabla operator, is the unit tensor and is the conformation tensor, which is symmetric and positive definite. The Reynolds number is described as Re = ( ℎ) ⁄ , where is the density, is the average flow velocity in the inflow and ℎ is the channel width. The Weissenberg number is defined as Wi = ℎ ⁄ , where is the relaxation time of the fluid. The viscosity ratio is defined as = + ⁄ = ⁄ , where and are the dynamic viscosity of the solvent and the solute of viscoelastic fluid, respectively.
is zero-shear rate viscosity. Unless otherwise specified, all data below are dimensionless.

Governing Equations and Numerical Method
Two-dimensional incompressible adiabatic flow was considered. The equations of conservations of mass and momentum were adopted. Generally, viscoelastic fluid has shear-thinning physical property, but to simplify the research complexity, Boger-type fluid with constant shear viscosity was selected in this paper, and the constitutive equation was Oldroyd-B [53]. This simplification is quite common for the convenience of research. [54,55]. For the sake of generality, each variable is nondimensionalized, and the dimensionless governing equations are the following.
where t is the time, u is the velocity, p is the pressure, ∇ is the nabla operator, I is the unit tensor and C is the conformation tensor, which is symmetric and positive definite. The Reynolds number is described as Re = (ρUh)/µ 0 , where ρ is the density, U is the average flow velocity in the inflow and h is the channel width. The Weissenberg number is defined as Wi = λU/h, where λ is the relaxation time of the fluid. The viscosity ratio is defined as β = µ s / µ s + µ p = µ s /µ 0 , where µ s and µ p are the dynamic viscosity of the solvent and the solute of viscoelastic fluid, respectively. µ 0 is zero-shear rate viscosity. Unless otherwise specified, all data below are dimensionless. As for the boundary conditions, the no-slip condition was imposed at walls, and the opening condition was applied at outlets. A fully developed velocity profile was set at the inlet to eliminate the influence of entrance. The elastic stress was relaxed at the initial moment, i.e., the conformation tensor was a unit tensor. The detailed description of numerical scheme and solution method are given in [56].

Grid Independence Validation
The standard T-shape geometry was used to verify the grid independence, and four sets of grids were adopted. Different numbers of uniform grids were arranged in the intersection center of the channel, and the same number of grid nodes were distributed in the three arms. Starting from the minimum mesh size, the mesh was distributed at a grid growth rate of 1.1 from the center to the outside. See Table 4 for the specific mesh settings. In this paper, the ranges of Wi and β were 0~10 and 0.1~0.9, respectively. An example with Wi = 10 and β = 0.1 was selected for grid independence verification. The time step for different mesh settings was 10 −4 and Re was 0.01. Note: ND x and ND y are the number of nodes in the x and y directions, respectively. ND branch is the number of nodes along the three branches. ∆x min or ∆y min is the minimum mesh size in the x or y direction. NC and ND are the total number of cells and nodes, respectively. Figure 2 shows that the statistical time-averaged flow rates of the two outlets tended to be stable with the increase of grid number, and the difference between the results of Mesh3 and Mesh4 was 1.2%. Considering the computing capacity and accuracy, the node layout of Mesh3 was selected during the subsequent research. It is worth noting that the results of the two branch flows of Mesh4 were opposite to those of Mesh3. For the case of Mesh4, Figure 3 shows that the flow rate at the right outlet was so large that the fluid flow in the left arm reversed to the right side, which was contrary to the results of the other three sets of grids. This random phenomenon in the flow direction was reasonable and consistent with the experimental [57] and numerical simulation [58] results of viscoelastic fluid under the microfluidic cross-slot channel, indicating the randomness of viscoelastic fluid in elastic instability. y direction. NC and ND are the total number of cells and nodes, respectively. Figure 2 shows that the statistical time-averaged flow rates of the two outlets tended to be stable with the increase of grid number, and the difference between the results of Mesh3 and Mesh4 was 1.2%. Considering the computing capacity and accuracy, the node layout of Mesh3 was selected during the subsequent research. It is worth noting that the results of the two branch flows of Mesh4 were opposite to those of Mesh3. For the case of Mesh4, Figure 3 shows that the flow rate at the right outlet was so large that the fluid flow in the left arm reversed to the right side, which was contrary to the results of the other three sets of grids. This random phenomenon in the flow direction was reasonable and consistent with the experimental [57] and numerical simulation [58] results of viscoelastic fluid under the microfluidic cross-slot channel, indicating the randomness of viscoelastic fluid in elastic instability. one exception of the last point corresponding to at the right outlet. The cross of the lines results from the random phenomenon in such flows as discussed above. To clarify this, the last two data are connected by dotted line.

Results and Discussion
Considering the influence of the initial conditions, the result analyses were performed for only the last 700 out of the whole 1000 dimensionless simulated time units.

Results and Discussion
Considering the influence of the initial conditions, the result analyses were performed for only the last 700 out of the whole 1000 dimensionless simulated time units. The pulsing performance of the outlet flow was examined, which was mainly measured by the mean and variance of the average flow. Then, the stretching degree of microstructure (represented by the trace of conformation tensor) and the flow field were studied to analyze the phenomenon. The Re for all cases was fixed to 0.01 but with various Wi, indicating that the flow phenomenon appearing at such conditions could be only ascribed to the elastic effect of viscoelastic fluid.

Pulsing Performance of Viscoelastic Fluid Flow in Standard T-Shaped (ST) Channel
Variables oscillating with time are the typical features for an oscillator. Therefore, we first discuss the change of the pressure and velocity. To study the flow characteristics of viscoelastic fluid, three typical Wi of 1, 5 and 10 and three typical viscosity ratios β of 0.1, 0.5 and 0.9, respectively, were selected as the research objects. After analyzing the results, it was found that the flow under some working conditions had regular oscillation characteristics. This subsection first shows the results of the case at Wi = 5 and β = 0.5 with expected results in the ST channel.
The variation of streamline velocity u* normalized by the maximum velocity at inlet and pressure at monitored points of (x, y) = (−5, 0) and (x, y) = (−10, 0), respectively, are shown in Figure 4. It can be clearly seen that for Newtonian fluid, the velocity and pressure at the monitoring point did not change with time, because the Re was so small (Re = 0.01) that the flow was in laminar state. For viscoelastic fluid with enough elastic effect, however, the flow velocity and pressure showed regular periodic oscillation. In addition, the fluctuation form of velocity developed from irregular periodic arrangement at the upstream to nearly perfect sinusoidal-like waveform at the downstream, indicating that the flow gradually fell into a fully developed state. This was also concluded from the velocity profile across the channel as drawn in Figure 5. The velocity profile of Newtonian fluid in laminar flow is a standard parabolic shape, while the velocity distribution of viscoelastic fluid in the upstream undeveloped state is deviated from the standard parabolic shape, and the location of the maximum velocity is 10% of the channel width away from the center toward the lower wall. The fully developed velocity profile of viscoelastic fluid at the downstream is almost coincident with that of Newtonian fluid, but still inclines slightly to the lower wall. the fluctuation form of velocity developed from irregular periodic arrangement at the upstream to nearly perfect sinusoidal-like waveform at the downstream, indicating that the flow gradually fell into a fully developed state. This was also concluded from the velocity profile across the channel as drawn in Figure 5. The velocity profile of Newtonian fluid in laminar flow is a standard parabolic shape, while the velocity distribution of viscoelastic fluid in the upstream undeveloped state is deviated from the standard parabolic shape, and the location of the maximum velocity is 10% of the channel width away from the center toward the lower wall. The fully developed velocity profile of viscoelastic fluid at the downstream is almost coincident with that of Newtonian fluid, but still inclines slightly to the lower wall.    The amplitude spectrum of axial velocity fluctuation of viscoelastic fluid in the ST structure is shown in Figure 6. There existed a fundamental frequency ( = 0.0571, the same value for the modified structures of CU1 and CU2) and two harmonic frequencies at (x, y) = (−5, 0). As the flow developed downstream, the amplitude of fundamental frequency increased, and that of the harmonic frequency weakened and some harmonic frequencies disappeared. The amplitude spectrum of axial velocity fluctuation of viscoelastic fluid in the ST structure is shown in Figure 6. There existed a fundamental frequency ( f 0 = 0.0571, the same value for the modified structures of CU1 and CU2) and two harmonic frequencies at (x, y) = (−5, 0). As the flow developed downstream, the amplitude of fundamental frequency increased, and that of the harmonic frequency weakened and some harmonic frequencies disappeared.

Effect of the Weissenberg Number and Viscosity Ratio
Given a constant average velocity unity at the inlet, the sum of the instantaneous average flow rate of the two outlets should be unity, considering the conservation of continuity. For simplicity, only the result of the left branch of the channel is considered in the following discussion. When the flow pattern is completely symmetric, the is equal

Effect of the Weissenberg Number and Viscosity Ratio
Given a constant average velocity unity at the inlet, the sum of the instantaneous average flow rate Q t of the two outlets should be unity, considering the conservation of continuity. For simplicity, only the result of the left branch of the channel is considered in the following discussion. When the flow pattern is completely symmetric, the Q t is equal to 0.5. Because of the possible complex flow, this parameter may exceed the range of 0 to 1, indicating that backflow may occur at an outlet.
The variation of instantaneous flow rate with time is shown in Figure 7.

Effect of the Channel Shape
The variation of the instantaneous flow rate at the left exit of the four modified structures with various symmetrical square grooves around the intersection of the T-shaped channel is shown in Figure 8. The general results were similar to those of the standard Tshaped channel as discussed in Section 3.2. However, the cavities carved on different positions exerted an obviously unique impact on the flow field. In the case of beta0.1-Wi10 in the green line, for example, the wider upper cavities (CU2 and CUD) seemed to have a suppression influence on the pulsing flow in terms of the amplitude. In addition, the results of beta0.5-Wi5, an example of medium viscosity with medium Wi, deserve special attention. Firstly, as shown in the subgraphs in Figure 8 and more clearly demonstrated

Effect of the Channel Shape
The variation of the instantaneous flow rate at the left exit of the four modified structures with various symmetrical square grooves around the intersection of the T-shaped channel is shown in Figure 8. The general results were similar to those of the standard T-shaped channel as discussed in Section 3.2. However, the cavities carved on different positions exerted an obviously unique impact on the flow field. In the case of beta0.1-Wi10 in the green line, for example, the wider upper cavities (CU2 and CUD) seemed to have a suppression influence on the pulsing flow in terms of the amplitude. In addition, the results of beta0.5-Wi5, an example of medium viscosity with medium Wi, deserve special attention. Firstly, as shown in the subgraphs in Figure 8 and more clearly demonstrated in Figure 9, for the standard T-type and the two cases of CU1 and CU2, the instantaneous average flow rate at the outlet presented approximately periodic fluctuation. However, for the two cases of CD and CUD with grooves at two lower corners, the instantaneous average flow rate at the outlet presented a weak fluctuation with a very small amplitude compared with the first three cases, as displayed in Table 5. Secondly, the closer the average value of the outlet flow rate approached to 0.5, the more likely the exits on both sides could be used as the pulsation excitation sources. In other words, the fan-out of the device could become two. In brief, considering both the pulsation intensity and the fan-out, the ST structure was the most suitable for the viscoelastic fluid microfluidic oscillator when β = 0.5 and Wi = 5 in this study. In addition, when the viscosity ratio was medium but the Wi was high, as shown with diamond symbols in Figure 8a, for example, the flow presented a certain regular pulsation, but the flow was easy to flip with different flow rates at the two outlets, as shown with diamond symbols in Figure 8b, for example, which was not easy to control.      To understand the motion of viscoelastic fluid under the conditions of β = 0.5 and Wi = 5, the variation of microstructure stretching in viscoelastic fluid flow was analyzed, which was tracked by the trace of conformation tensor tr(C) [59]. Figure 10 shows the distributions of tr(C) for various geometries at the instant with the largest flow rate in Figure 9. In the ST channel, the relatively large tr(C) region dominated by shear stress tended to be close to the vertical axis of symmetry at the throat in the inflow channel, while after the flow entered the junction, the macromolecular structure was stretched greatly at the stagnation point dominated by tensile stress, and a severe stretching zone was formed near the upper wall. The situations of cases CU1 and CU2 were similar to that of ST. However, for the geometries with square recesses added to the lower walls (CD and CUD), the relatively large tr(C) region tended to be far away from the vertical axis of symmetry at the throat due to the smaller pairs of vortices, as shown in Figure 12. In addition, compared with the structures without cavities at the two lower corners (ST, CU1 and CU2), the range of low tr(C) region increased in the junction due to the presence of sudden expansion structure. Figure 11 shows the variation trends of time-averaged tr(C) along the two axes. It can be seen from Figure 10 that the position where the maximum tension occurred in the horizontal channel was close to the upper wall, which had an important influence on the results. Therefore, Figure 11a shows the variation trend of tr(C) on the horizontal line with y = 0.4. All the curves are symmetrical about the vertical axis, and only the left halves are plotted. The stagnation points of ST and CD structures are fixed on the intersection of the horizontal wall and the vertical axis, and their tr(C) takes the maximum values locally. However, for the structures of CU2 and CUD with grooves on the upper side, their stagnation points are free to move into the groove, as shown in Figure 10c,e. Therefore, the maximum values of tr(C) are located in the groove, and the corresponding curves on the axis of symmetry in Figure 11a take the local minimum value. It is worth noting that the width of the groove also had an important effect on the flow. The results of the CU1 structure with smaller groove width were similar to those of the ST structure, which seemingly indicated that the smaller groove size had almost no effect on the flow. The variation of tr(C) on the centerline of axis y shown in Figure 11b shows again that the maximum tensile position of the structure with a wider groove on the upper side is inclined to the side of the upper groove.
distributions of tr(C) for various geometries at the instant with the largest flow Figure 9. In the ST channel, the relatively large tr(C) region dominated by she tended to be close to the vertical axis of symmetry at the throat in the inflow while after the flow entered the junction, the macromolecular structure was greatly at the stagnation point dominated by tensile stress, and a severe stretch was formed near the upper wall. The situations of cases CU1 and CU2 were simil of ST. However, for the geometries with square recesses added to the lower walls CUD), the relatively large tr(C) region tended to be far away from the vertical axi metry at the throat due to the smaller pairs of vortices, as shown in Figure 12. In compared with the structures without cavities at the two lower corners (ST, CU2), the range of low tr(C) region increased in the junction due to the presence o expansion structure. Figure 10. Distribution of the trace of conformation tensor tr(C) for various geometries at the instant with the largest flo Figure 9. For all cases, = 0.5, Wi = 5. Figure 11 shows the variation trends of time-averaged tr(C) along the two ax be seen from Figure 10 that the position where the maximum tension occurred in izontal channel was close to the upper wall, which had an important influenc results. Therefore, Figure 11a shows the variation trend of tr(C) on the horizontal y = 0.4. All the curves are symmetrical about the vertical axis, and only the left h However, for the structures of CU2 and CUD with grooves on the upper side, their stagnation points are free to move into the groove, as shown in Figure 10c,e. Therefore, the maximum values of tr(C) are located in the groove, and the corresponding curves on the axis of symmetry in Figure 11a take the local minimum value. It is worth noting that the width of the groove also had an important effect on the flow. The results of the CU1 structure with smaller groove width were similar to those of the ST structure, which seemingly indicated that the smaller groove size had almost no effect on the flow. The variation of tr(C) on the centerline of axis y shown in Figure 11b shows again that the maximum tensile position of the structure with a wider groove on the upper side is inclined to the side of the upper groove.

Mechanism Analysis
As described above, the viscoelastic flows in the ST channel and channels with symmetrical modifications on the upper wall (CU1 and CU2) exhibited approximately sinusoidal-like oscillation behaviors at the two outlets. The local flow conditions near the junction were responsible for the oscillatory flow downstream the exit channels. As depicted in Figure 12, for the geometries without cavities at the lower corners (ST, CU1 and CU2),

Mechanism Analysis
As described above, the viscoelastic flows in the ST channel and channels with symmetrical modifications on the upper wall (CU1 and CU2) exhibited approximately sinusoidallike oscillation behaviors at the two outlets. The local flow conditions near the junction were responsible for the oscillatory flow downstream the exit channels. As depicted in Figure 12, for the geometries without cavities at the lower corners (ST, CU1 and CU2), a pair of lip vortices formed in the entrance branch at the throat region moved up and down periodically just as the periodically changing flow rate at outlets shown in Figure 9. For the cases of CD and CUD, only a pair of small vortices were generated in the throat region, which were not strong enough to impart significant effect on the flow field in the downstream, resulting in a relatively steady flow both in the mainstream and in the cavities. Also, it can be seen from Figure 12 that the velocity magnitude at the lower part of the horizontal branches was higher than that at the upper part, which directly explained the fact that the undeveloped velocity profile as shown in Figure 5 bent to the lower wall. grooves on the upper side (CU1, CU2 and CUD). Here, we used L (the distance between the second and third maxima of tr(C), as shown in Figure 14) to measure the elastic effect in the junction domain. The smaller L showed a stronger elastic effect in the whole junction domain. As shown in Figures 13 and 14, it indicated that a T-shaped-based structure with a smaller L in the junction tends to be more efficient as a microfluidic oscillator.  The viscoelastic fluid in the entrance branch was only subjected to strong shear action near the wall. As illustrated previously, the dynamics of the vortex pair in the throat region played an important role. Therefore, we took the variation of tr(C) along the line x = −0.4 (close to the wall) to analyze the microstructure stretching quantitatively as displayed in Figure 13. The stretching of microstructures induced the peculiar flow behaviors. It can be seen from Figure 13 that, along this line close to the wall, the stretching extent of the molecular structure increased and decreased three times alternately. The first minimum was related to the vortex at the upstream of the sudden expansion, while the second was due to the weak shear effect of the main flow near the center of the horizontal channel, and the last was related to the low-speed motion near the upper wall (ST and CD) or the grooves on the upper side (CU1, CU2 and CUD). Here, we used L (the distance between the second and third maxima of tr(C), as shown in Figure 14) to measure the elastic effect in the junction domain. The smaller L showed a stronger elastic effect in the whole junction domain. As shown in Figures 13 and 14, it indicated that a T-shaped-based structure with a smaller L in the junction tends to be more efficient as a microfluidic oscillator. To further verify that the oscillation characteristics were related to the periodic fluctuation of the upstream vortex pair and the distance L, viscoelastic fluid flow in a channel with a groove carved inward on the upper wall (denoted as case CU3, inspired by the reference [60]) was also simulated under the same condition with = 0.5 and Wi = 5, as shown in Figure 15 for typical results. Half-channel width was chosen as the characteristic length for CU3. The results showed that the standard deviation of the outlet flow To further verify that the oscillation characteristics were related to the periodic fluctuation of the upstream vortex pair and the distance L, viscoelastic fluid flow in a channel with a groove carved inward on the upper wall (denoted as case CU3, inspired by the reference [60]) was also simulated under the same condition with = 0.5 and Wi = 5, as shown in Figure 15 for typical results. Half-channel width was chosen as the character- To further verify that the oscillation characteristics were related to the periodic fluctuation of the upstream vortex pair and the distance L, viscoelastic fluid flow in a channel with a groove carved inward on the upper wall (denoted as case CU3, inspired by the reference [60]) was also simulated under the same condition with β = 0.5 and Wi = 5, as shown in Figure 15 for typical results. Half-channel width was chosen as the characteristic length for CU3. The results showed that the standard deviation of the outlet flow rate σ(Q t ) was increased by 129% (calculated by (σ CU3 and σ(Q t ) ST are the standard deviation of Q t for cases of CU3 and ST, respectively) from 0.0075 to 0.0172, compared with the case ST. Note that the averaged flow rate Q t at the left outlet was 0.4661, deviated from the ideal value of 0.5. As drawn in Figure 14, the distance L of case CU3 is dramatically decreased compared to that of the case ST. It is worth noticing that the oscillation frequency of the case CU3 increased sharply from 0.0571 of the case ST to 0.1250, implying that the oscillating performance is closely related to the channel geometry in addition to the physical properties of viscoelastic fluid.  The results of viscosity ratio β = 0.9 were similar to those of Newtonian fluid, therefore the corresponding results are not presented in this paper. In sum, the pulsation intensity of the ST and its modified geometry in this study were at most only 3.7% of the mean value. For now, regarding the practical application of the so-called microfluidic oscillator as an "excitation source" in microfluidic circuits, there seems to be a certain distance between the research and reality. The geometry form of CU3 that had an obvious improved performance was not optimized systematically, but was a preliminary structure obtained by mechanism analysis. Inspired by the elegant and ingenious structure proposed by predecessors (for example, a variant of switch structure in [28]), we can expect to get a better performance of the simple T-shaped microfluidic oscillator with viscoelastic fluid by carefully designing the core components, together with three-dimensional simulation and later experimental verification.

Conclusions and Outlook
In this paper, the flow of Boger-type viscoelastic fluid in two-dimensional standard T-shaped and its modified structures was numerically simulated. The main conclusions are drawn as follows: 1. For viscoelastic fluid with medium viscosity ratio and medium Wi the average flow rate at the outlets of standard T-shaped and its modified structures changes approximately sinusoidal-like with time; 2. To generate regular periodic signals, both Wi and viscosity ratio need to be within a Figure 15. The variation of instantaneous flow rate Q t at the left outlet with dimensionless time for the case of standard T-shaped channel (ST) and modified channel with cavity on the top wall with half-channel width carved out (CU3). For all cases, β = 0.5, Wi = 5. Inset: Three instants for contour plots of velocity magnitude u superimposed with streamlines. T is the period corresponding to the fundamental frequency of the CU3 structure.
The results of viscosity ratio β = 0.9 were similar to those of Newtonian fluid, therefore the corresponding results are not presented in this paper. In sum, the pulsation intensity of the ST and its modified geometry in this study were at most only 3.7% of the mean value. For now, regarding the practical application of the so-called microfluidic oscillator as an "excitation source" in microfluidic circuits, there seems to be a certain distance between the research and reality. The geometry form of CU3 that had an obvious improved performance was not optimized systematically, but was a preliminary structure obtained by mechanism analysis. Inspired by the elegant and ingenious structure proposed by predecessors (for example, a variant of switch structure in [28]), we can expect to get a better performance of the simple T-shaped microfluidic oscillator with viscoelastic fluid by carefully designing the core components, together with three-dimensional simulation and later experimental verification.

Conclusions and Outlook
In this paper, the flow of Boger-type viscoelastic fluid in two-dimensional standard T-shaped and its modified structures was numerically simulated. The main conclusions are drawn as follows:

1.
For viscoelastic fluid with medium viscosity ratio and medium Wi the average flow rate at the outlets of standard T-shaped and its modified structures changes approximately sinusoidal-like with time; 2.
To generate regular periodic signals, both Wi and viscosity ratio need to be within a certain range, beyond which no oscillation or completely irregular flow will occur; 3.
The mechanism of the oscillation characteristics is related to the periodic fluctuation of the upstream vortex pair and the size of the elasticity-dominated area in the whole junction domain. With the presence of the dynamic evolution of the pair of vortices in the upstream near the intersection, the oscillating intensity increases as the elasticitydominated area in the junction enlarges; 4.
When the T-type simple structure is considered as the potential realization of the oscillator, for now, the modified structure with a groove carved inwards on the upper wall facing the entrance branch is the most suitable for the oscillator to provide excitation for the downstream equipment.
At present, the results obtained above are based on the constitutive equation of the Oldroyd-B model, which is suitable for the case where the stretching length of molecules is much less than the maximum extensibility. The extent of microstructure stretching (measured by the trace of conformation tensor tr(C)) for the considered ranges of Wi and β in various geometries in this paper reached values in the order of 2000, as displayed in the Figure 10. Given the more accurate simulation when the stretching length tends to the maximum stretching length under consideration, the application of the FENE model to carry out related research work will also be necessary in the future. Moreover, only sinusoidal-like periodic pulsating flow with limited amplitude is generated with two-dimensional simulation. To make full use of the T-type channel, which has the advantages of simplicity, no moving parts and fan-out of two, a lot of work still needs to be done, including ingenious design of the channel shape and careful selection of the appropriate flow condition range. As a next step of research, a thorough three-dimensional simulation using a realistic fluid model with finite extensibility such as the FENE model and meticulous microchannel experiments using real viscoelastic fluids are required to realize a T-type microfluidic oscillator with considerable pulsation amplitude to excite the downstream flows.