1 Introduction

In a fluid domain,flow parameters can vary with time and space. These variations can lead to streamline bending and pulsation,resulting in vortex formation and turbulence as shown in Fig. 1^[1]. Reynolds’ experiments^[2] revealed the transformation of flow regime from laminar to turbulent flow. Ludwig Plandtl^[3] proposed the boundary-layer theory. Boundary layer is defined as a thin layer close to a wall surface where the velocity gradient is large. Flow is considered to be potential outside the boundary layer. The boundary layer theory defines the transition of flow regime from turbulent flow in the boundary layer to laminar flow outside the boundary layer. In this study,the emergence of vortices in a gradient flow distribution is investigated. In the previous studies,several theoretical investigations and experimental measurements have been conducted to define the vortex formation close to a wall surface (Mc-Quivey and Richardson^[4],Esfahanian et al.^[5],Li et al.^[6],Summer and Akosile^[7],Wang et al.^[8],Huang and Zhou^[9],Ehrenstein et al.^[10],Gires et al.^[11],K´ad´ar and Balan^[12],and Ashrafi and Hazbavi^[13]).

In this study,for the first time,differential equations are solved considering the inertia term along the flow direction. The Oseen transformation is used to decompose the differential equations,and the vorticity,velocity,and pressure solutions are obtained using a new contour integral method. Theoretical results can be used in a variety of problems related to vortex formation close to the wall surface. 2 Vorticity equation of linear gradient flow distribution

A linear gradient flow distribution close to a flat wall is studied as shown in Fig. 2. The flow direction is along the x-direction. The maximum velocity at the upper boundary is given as U. The velocity is expressed as u = z = kz,where k = is the velocity gradient,and h is the total depth. The time-averaged flow velocity is given as V = (kz,0,0). The vortex flow velocity vector V is given as V = (u,v,w). The instantaneous flow velocity vector V ′ is given as V ′ = (u′,v′,w′),where V ′ = V + V . The vorticity in the y-direction is defined as w_y = . The instantaneous flow satisfies the incompressible viscous equation of continuity and momentum. In the Cartesian coordinates,these equations are expressed by

where t is the time, is the gradient operator,ρ is the fluid density,p′ is the instantaneous pressure,ν is the kinematic viscosity coefficient, is the Laplacian operator,and F is the body force. Using the Oseen transformation^{[14, 15]},the instantaneous velocity V ′ and the instantaneous pressure p′ can be represented as

Using Eq. (3) and neglecting the terms related to the second and higher orders of small parameters,Eqs. (1) and (2) can be rewritten as follows:

Looking at non-dimensional parameters,we can consider the Reynolds number as Re =,where Fr is the Froude number,the dimensionless velocity is ,and the dimensionless pressure is . The dimensionless coordinates are . The dimensionless time is . The dimensionless averaged velocity is . The corresponding equations of each εth-order are

Substituting the dimensionless parameters into Eq. (5),we have the hydrostatic pressure equation as

Substituting the dimensionless parameters into Eq. (6) and removing the brackets for simplicity,we have

Taking divergence of both sides of the second equation of Eq. (8) and using the continuity equation (the first equation of Eq. (8)),we have

Therefore,the ε¹th-order pressure function can satisfy the Poisson equation and is directly proportional to the vertical velocity rate along the flow direction. Taking the curl of both sides of the second equation of Eq. (8) and the continuity equation (the first equation of Eq. (8)),we have

Assuming a vertical two-dimensional problem,we have ν = 0 and = 0. Equation (10) can be simplified as

We consider an analytic function in the form of

Then,we have

The above equation is the ε1th-order vorticity function differential equation. 3 Solution to vorticity function equation

We assume that Eq. (13) has a variable separable solution in the form of

Substituting Eq. (14) into Eq. (13) and considering the quasi-periodicity of T (t),we have

where ω is an undetermined constant,and i is the imaginary unit. Equation (15) can be presented as

In Eq. (16),the solution to the first equation is

where A is a constant. The second equation has quasi-periodicity variations along the xdirection. The solution is assumed to be in the form of where n = is the wave number,σ = 1,2,3,· · · ,and L is the wavelength. We can obtain a set of solutions corresponding to each value of n. Without loss of generality,we can solve the equation for any value of n. Substituting Eq. (18) into the second equation of Eq. (16) yields

Let K_n = and ξ_n = −(n² + iωRe)K_n^-2 − iK_n^z . Then,Eq. (19) can be rewritten as

Equation (20) is the complex variable form of the Airy equation^[16]. Using ξ_n = −ξ′_n transform,we can obtain the equivalent equation as

4 Contour integral and solution to complex variable Airy equation 4.1 Contour integral

We consider a complex analytic function as

For the contour Γ integral,w and ξ_n are complex variables. The contour Γ starts from the origin O and goes to A along the real axis R. The length of OA is ρ (ρ > 0). The arc continues to B,which has the position of ρe^iπ/6. Then,it returns to the origin O along the straight-line BO. Then,it starts from the origin O and goes towards C along a straight line,which has the position of ρe^i5π/6. After that,the arc goes to D,which is on the imaginary axis at the position iR,as shown in Fig. 3. According to the Cauchy integral theorem along the closed curve^[17], f(w) can be resolved in the contour. Therefore,the integral is zero along the contour Γ,i.e.,

The following part studies the solution for each integral section. In order to compute the OA integral section,which is along the real axis,we consider w = τ,where τ is a real number. We have

In order to compute the integral section,which is along an arc of radius ρ,we consider w = ρe^iθ and dw = iρe^iθdθ. We have

The integral is independent of ρ. Then,we substitute e^i3θ = cos(3θ) + i sin(3θ) into the above equation (25),and we have

When ρ → ∞,the convergence of Eq. (25) is related to . When θ varies between 0 and ,3θ varies between 0 and . Therefore,cos(3θ) > 0. We have

Figure 4 shows that the amplitude of changes with ρ. When , and tends to 0 faster than . We compute the BC section integral,which is along the straight line of τe^iπ/6,by letting w = τe^iπ/6 = τe^{iπ(2/3−1/2)} = −iτe^2π/3. We have

We compute the OC section integral,which is along the straight line of τe^i5π/6,by letting w = τe^i5π/6 = −τe^iπ/6 = −iτe^−i2π/3. We have

We compute the integral,which is along an arc of radius ρ,by letting w = ρe^iθ and dw = iρe^iθdθ. We have When θ varies between and ,3θ varies between and ,and cos(3θ) ≥ 0. Therefore, we have

We compute the DO integral,which is along the imaginary axis i by letting w = τe^iπ/2 = iτ. We have From Eq. (24) to Eq. (28),when ρ → ∞,both and integrals tend to zero. Then,by looking at Eq. (23),we have an integral summation along the path ρe^iπ/6 ~ 0 and the path 0 ~ ρe^i5π/6. 4.2 Solution to complex variable Airy equation

We assume the solution to the complex variable Airy equation in the form of

where H_nj is an undetermined constant,and j = 1,2. We substitute Eqs. (31) and (32) into the complex variable Airy equation. Noticing that 2π/3 = π/2 + π/6,we have

Similarly,we can prove

The above equation (33) is the solution to the equivalent Airy equation (21). The results show that the left side of Eq. (30) satisfies the complex variable Airy equation. Therefore, we can prove that the right side of Eq. (30) is also the solution to the complex variable Airy equation. We can write the general solution form as

We assume a solution in the form of The general solution to the complex variable Airy equation can be written as

Substituting the solution into Eq. (14),we can obtain the solution of nth-order vortex function as

where H_n1,H_n2,n,and ω can be determined by the boundary conditions and the physical conditions. 5 Velocity solution of pulsatile flow field 5.1 Velocity equation

According to the vorticity equation and the continuity equation,we have

Then,we can have Considering = −iK_n and Eq. (37),we have where 5.2 Homogeneous solution and special solution of horizontal velocity component

We assume u = U(z) exp(i(nx + ωt)) and substitute it into Eq. (40). We have

The homogeneous solution to Eq. (45) is in the form of

where D_n1 and D_n2 are constant,which are determined by the boundary conditions. We assume the special solution to Eq. (45) in the form of where G_x(τ),K_x(τ),Q_x(τ),and S_x(τ) are complex functions. Substituting U₂ into Eq. (45) yields If we compare the coefficients,we have We can obtain the special solution of horizontal velocity as where the τ integral is from 0 to ∞,and there exists a singularity. Therefore,we need to solve this problem by means of singularity contour integral. 5.3 Singularity contour integral of horizontal velocity

We assume a complex function in the form of

In the τ plane,the contour starts from the origin O,goes to A,and then goes to C passing B. ABC is a semicircle arc of radius ε. Here,we bypass the singularity . Then,we go from C to D along the τ-axis. From D,an arc of radius ρ in the counter-clockwise direction reaches the point E. Then,from E,we go back to O along . Then,we go from the origin O to F along with an arc of radius ρ in the clockwise direction of the point G on the iτ-axis. Finally,we go from G to the origin O along the iτ-axis,as shown in Fig. 5. Along the contour ,the integral is zero,and we have

In the straight part O → E,we have w = τ and dw = dτ.

In the curved arc ,we have w = + εe^iθ and dw = iεe^iθdθ. We have When ε → 0,we use L’Hospital’s rule and we have

In the curved arc ,we have w = ρe^iθ and dw = iρe^iθdθ.

where . When ρ → ∞,the convergence of Eq. (55) depends on . When θ varies between 0 and ,3θ varies between 0 and π, and sin (3θ) ≥ 0. Therefore,we have

For the straight part EO,we have w = τe^iπ/3,dw = e^iπ/3dτ,and Using the angular transformation,we have

For the straight part OF,we have w = τe^-iπ/3,dw = e^-iπ/3dτ,and We use the angular transformation and have

For the curved arc ,we have w = ρe^iθ,dw = iρe^iθdθ,and where . When ρ → ∞,the convergence of Eq. (57) depends on . When θ varies between ,3θ varies between −π and ,and sin(3θ) ≥ 0. Therefore,we have

For the straight part GO,we have w = −iτ,dw = −idτ,and Then,we substitute Eqs. (53)-(59) into Eq. (52),and we solve the integral terms of Eq. (50), For simplicity,we use ,where i and j are perpendicular unit vectors. L_x(ax) can be extended to any dimensional vector. In this study, a = e^i2π/3,b = (K_nτ)²e^−i2π/3− n²,c = (K_nτ)²e^i2π/3− n²,and x = −i(±ξ_nτ). The properties of L_x(ax) derivatives are as follows:

Circulation will appear after solving the third derivative since a² = e^i4π/3 = e^{i(2π−2π/3)} = e^−i2π/3 = a⁻¹ and L_x′′′(ax) == L_x(ax). We substitute Eq. (60) into Eq. (50), and the velocity special solution in the horizontal direction is obtained as

The velocity general solution in the horizontal direction is

where N = −(n² + iωR)K_n² and M = n/K_n. 5.4 General solution and special solution of vertical velocity component

We assume w = W(z) exp(i(nx + ωt)) and substitute w into Eq. (40). We have

The homogeneous solution to Eq. (63) can be presented as where E_n1 and E_n2 are undetermined coefficients. The special solution to Eq. (63) can be presented as where G_z(τ),K_z(τ),Q_z(τ),and S_z(τ) are the complex functions. Then,substituting the special solution W₂ into Eq. (63),we have Comparing the coefficients,we have We can obtain the special solution of velocity in the vertical direction as The above equation (68) contains singularity. Therefore,we need to solve this problem by means of singularity contour integral. 5.5 Singularity contour integral of vertical velocity

We assume an analytic complex function as

The contour integral is computed term by term using the same method discussed in Section 5.3. Therefore,Eq. (68) can be presented as

The vertical velocity is presented as 6 Pressure function in shear flow field

In the Couette flow problem,the pressure gradient is used to find the velocity distribution. However,in the study,the velocity distribution is used to find the pressure. Rewrite Eq. (71) as

We consider the pressure in the form of p = P (z) exp(i(nx+ωt)) and substitute p into Eq. (9). We have

The homogeneous solution to Eq. (73) is presented as where F_n1 and F_n2 are undetermined coefficients. The special solution to Eq. (73) is presented as where G_n1 and G_n2 are undetermined coefficients. We consider L_xg(ax) = (g₁(τ)i − g₂(τ)j) · (e^axi+e^a−1 j),and L_xk(ax) = (k₁(τ)i −k₂(τ)j) · (e^axi+e^a−1xj). g₁(τ),g₂(τ),k₁(τ),and k₂(τ) are also undetermined functions. Substituting the special solution P₂ into Eq. (73),we have Comparing the coefficients,we have Considering both the homogeneous solution and the special solution,the following equation for the pressure function can be presented: where 7 Boundary conditions for finite problems

In order to determine the undetermined constants in Eqs. (37),(62),(71),and (81),the boundary conditions of the vorticity solution,the velocity solution,and the pressure solution are discussed in the following sections. The boundary conditions have various forms,but determining the boundary conditions of the complex problems has never been an easy task. A simple form of the boundary condition is presented in the following sections,and the variations of the solutions are discussed. 7.1 Vorticity solution

In the vorticity solution,the complex variable Airy functions A_j(ξ_n) and B_j(ξ_n) satisfy the complex variable Airy equation. When ξ_n is determined,the complex variable Airy function cannot converge at the same time,if ξ_n = −a − ib,a > 0,and b > 0,where both a and b have finite values. The integral ∫₀^∞ expdτ converges. However,the convergence of the integral i is determined before ξ_n. We consider positive sign as divergence and negative sign as convergence,respectively. Therefore,A_j(ξ_n) diverges,and B_j(ξ_n) converges. Removing the part related to A_j(ξ_n) and letting H_n1 = 0,we can obtain the vorticity solution as

The boundary condition φn|_z=0 = f₁ at the wall surface can determine H_n2. f₁ is the vorticity, which is induced by a rough wall surface,a velocity gradient,etc. In general,the above equation can be rewritten as From the above equation,it can be seen that after determining the time and the x-direction spatial position,B_j(ξ_n) can determine φ′_n variation. Using the previous studies,we can determine the parameter values. We obtain the parameter R,considering the turbulent boundary layer thickness formula 11.6^[18],where v_* is the friction velocity,δ is the boundary layer thickness,and ν is the kinematical viscosity coefficient. The circular frequency is considered to be equal to ω = 18.84 s⁻¹,which is obtained by the turbulent energy. The vortex wavelength is L = 2πh. The wave number is n = σ/h,where the relative depth is h = 1.0 m and σ = 1,2,3,· · · ,100. Therefore,we can obtain the variation of the complex variable parameter ξ_n as shown in Fig. 6. When the value of n increases,the real part of ξ_n increases,the imaginary part of ξ_n approaches a stable value,and the water depth of imapinary part decreases. Figures 7 and 8 show the variation of B_j(ξ_n) along the vertical direction. We can see that the value of B_j(ξ_n) is larger at the region close to the wall surface. This implies that when the value of n increases,the value of B_j(ξ_n) becomes larger.

7.2 Velocity solution

In order to find the velocity solution,we look at Eqs. (62) and (71) and have

The bottom velocity boundary conditions are u|_z=0 = 0 and w|_z=0 = 0. The upper velocity boundary conditions are u|_z=h = 0 and w|_z=h = 0. Then,we have and

We substitute the vorticity solution parameters into Eqs. (90) and (91),and we obtain the velocity distribution along the vertical direction (as shown in Figs. 9 and 10).

It can be observed from Figs. 9 and 10 that the horizontal velocity distribution along the vertical direction is consistent with the vertical velocity distribution. When n is small,the fluctuation of the velocity amplitude is vertically symmetric. We can observe a downward shift of the curve centroid with the increase in the value of n. In this case,the upper curvature decreases and the bottom curvature increases when n > 7. When n > 10,the bottom of the curves appears to be folded. That is the reason for the vortex formation close to the wall. 7.3 Pressure solution

Substituting the constant values determined by the vorticity solution into Eq. (81) yields

The upper boundary condition for the relative pressure is p|_z=h = 0. We assume that the bottom relative pressure is zero without the vorticity disturbance,which gives p|_(z=0,Hn2=0) = 0. Then,we have

and If we know the time and the x-direction spatial position and substitute the vorticity solution parameters into Eq. (95),we can obtain the pressure function distribution along the vertical direction (as shown in Fig. 11). When n is small,the pressure presents a linear distribution. When n = 9,the upper part of curve is mostly convex. When n > 9,the curve swings toward the negative direction. When n = 11,the negative pressure amplitude appears near the bottom. 8 Stability analysis 8.1 Vorticity distribution

The vorticity spatial distribution can be obtained using Eq. (83) as shown in Fig. 12. It can be observed from Fig. 12 that a larger value of vorticity appears at a region close to the bottom of the wall. The vorticity has an alternatively positive-negative pattern,when the relative water depth is h ≥ 0.3 and the relative vorticity is |φ′| < 0.1.

8.2 Turbulent motion stability

The particle relative trajectory formula can be obtained from Eqs. (90) and (91),i.e.,

We consider the initial position as (x₀,z₀)=(0.0,0.2). Here,we study two cases,i.e.,n=1 m⁻¹ and n=10 m⁻¹. Figures 13 and 14 show the particle relative trajectories for low wave numbers and high wave numbers,respectively. Figures 13(a)-13(d) show the mean flow equal to U_0.2 = 0.0,U_0.2 = 0.5u_max,U_0.2 = 1.0u_max,and U_0.2 = 3.0u_max,respectively,where U_0.2 denotes the mean velocity at the relative water depth equal to 0.2,and umax denotes the maximum horizontal velocity of particle. It can be observed from Fig. 13 that when the mean velocity exists,the particle trajectory is more symmetric. When no mean velocity exists,the particle trajectory takes the shape of an ellipse. The waveform shows symmetrical distribution along the x-direction and has stable shape.

It can be observed from Fig. 14 that,when the wave number is large,the particle trajectory has both asymmetric shape and inclination. When no mean velocity exists,the particle trajectory is inclined in the form of an ellipse. When the mean velocity increases,the waveform unfolds successively. It can be observed that larger wave number leads to earlier stability loss. 8.3 Distribution of pressure function

Figure 15 shows the pressure function (see Eq. (95)) distribution in a vertical plane. It can be observed that a larger value of the relative pressure appears at the region close to the bottom of the wall. The relative pressure has an alternatively positive-negative pattern closer to the bottom of the wall. Changes of the positive-negative amplitude are more violent when the relative water depth is h ≥ 0.4 and the relative pressure is |p′| < 0.1.

9 Comparison between theoretical results and experimental data 9.1 Comparison of velocity data

The velocity fluctuations can be presented as different frequency components using statistical average methods. In this paper,the velocity fluctuation is presented using Eqs. (90) and (91) and a simple combination solution,

where A_i and B_i (i = 1,2,3,4) are undetermined parameters,and n₁ and n₂ are two sets of different wave numbers. Then,we have

Using the boundary conditions,we have For instance,in the u-direction,we have In order to make Eq. (101) have nonzero solution,we should have λ = −1. Therefore,we have . If we consider a=-A₄,we have Similarly, The distribution curves of can be obtained by selecting the appropriate values of parameters a,n₁,n₂ and a′,n′ ₁,n′₂. Figures 16 and 17 show the comparison between the theoretical results of the velocity fluctuation and the measured data. The measured data are obtained from Prashun,Raichlen,Rao,and McQuivey’s experiments results^[19],which is Case 1. Figures 18 and 19 show the comparison between the theoretical results of the velocity fluctuation and the measured data. The measured data are obtained from Grass’ s experiments results^[19],which is Case 2. Table 1 shows the selected parameters.



9.2 Comparison of theoretical results and experimental data of turbulent shear force

We consider the turbulent shear force formula as

The following equation is obtained by the same process as Eqs. (102) and (103):

Figures 20 and 21 are prepared using the parameters provided in Table 2. It can be observed from Fig. 16 to Fig. 21 that the variation of experimental data agrees well with that of the theoretical results.

10 Conclusions

(i) This paper uses the Navier-Stokes equations and applies the Oseen transformation in order to find the velocity fluctuation in a finite problem using the vorticity function,the velocity function,and the pressure function in a shear flow field.

(ii) The complex solution of the Airy function is studied using a new contour integral method. Then,the vorticity function,the velocity function,and the pressure function are analytically solved in a vertical linear gradient flow distribution.

(iii) This paper presents an L_x-function form of the third derivative circulation to simplity the solution.

(iv) The vorticity and pressure in a linear gradient flow distribution are computed for given boundary conditions. The stability of particle motion trajectories is analyzed.

(v) A comparison between the theoretical results and the experimental measurements is provided with satisfactory agreement.

Acknowledgements The authors thank the State Key Laboratory of Hydraulic Engineering Simulation and Safety (Tianjin University) for support and thank Prof. R. A. DALRYMPLE of Johns Hopkins University for assistance.