Design of viscoelastic coatings to reduce turbulent friction drag

Abstract

A method is provided to select appropriate material properties for turbulent friction drag reduction, given a specific body configuration and freestream velocity. The method is based on a mathematical description of the balance of energy at the interface between the viscoelastic surface and the moving fluid, and permits determination of the interaction of turbulent boundary layer fluctuations with a viscoelastic layer by solving two subtasks—i.e., a hydrodynamic problem and an elasticity problem, which are coupled by absorption and compliancy coefficients. Displacement, velocity, and energy transfer boundary conditions on a viscoelastic surface are determined, and a Reynolds stress type turbulence model is modified to account for redistribution of turbulent energy in the near-wall region of the boundary layer. The invention permits drag reduction by a coating with specified density, thickness, and complex shear modulus to be predicted for a given body geometry and freestream velocity. For practical applications, viscoelastic coatings may be combined with additional structure, including underlying wedges to minimize edge effects for coatings of finite length, and surface riblets, for stabilization of longitudinal vortices.

Description

parts of the solution are coupled by boundary, conditions,
An overbar indicates time-averaging:
U_i=Ū_i  (Equations 8a-8b-8c)
Substituting Equations 7a-7c into Equations 6a-6c and time-averaging yields the following system of three complex nonlinear second-order partial differential equations of motion for turbulent flow: $\begin{array}{cc}{U}_j\frac{\partial U_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial 𝒫}{\partial x_i}+g_i+v\frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}-\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_j}& \left(\mathrm{Equations}9a\text{-}9b\text{-}9c\right)\end{array}$
In Equations 9a-9c, the components u′², v′², and w′² are termed normal Reynolds stresses, and the components in the form − u′v′, − v′w′, and − u′w′ are termed Reynolds shear stresses. The turbulent kinetic energy, k, is defined as: $\begin{array}{cc}k=\frac{1}{2}\left(\stackrel{\_}{u^{\mathrm{\prime 2}}}+\stackrel{\_}{v^{\mathrm{\prime 2}}}+\stackrel{\_}{w^{\mathrm{\prime 2}}}\right)& \left(\mathrm{Equation}10\right)\end{array}$
Closure of the generalized system of equations including the continuity equation (Equation 5) and the equations of motion (Equations 9a-9c) for a turbulent flow requires seven additional equations to characterize the six Reynolds stresses, u_i′u_j′, and the rate of transport of the turbulent kinetic energy, k. This invention solves for the isotropic dissipation rate, ε, which is related to energy transport through the fluid and at the fluid-surface interface. The energy transport equation is based upon the first law of thermodynamics, where heat, dQ, added to a volume during an element of time, dt, serves to increase internal energy, dE, and to perform work, dWK. $\begin{array}{cc}\frac{dQ}{\underset{\mathrm{heat}}{\underbrace{dt}}}=\frac{dE}{\underset{\mathrm{energy}}{\underbrace{dt}}}+\frac{d\mathrm{WK}}{\underset{\mathrm{work}}{\underbrace{dt}}}& \left(\mathrm{Equation}11\right)\end{array}$

There exist multiple approaches within the literature for developing additional equations for Reynolds stress terms in turbulent flow, but this invention adopts a Reynolds-stress-transport-type methodology. In this methodology, equations for Reynolds stresses take the following general form: $\begin{array}{cc}{U}_k\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}=P_{\mathrm{ij}}-{\Pi }_{\mathrm{ij}}-\frac{\partial J_{\mathrm{ijk}}}{\partial x_k}-2{\varepsilon }_{\mathrm{ij}}& \left(\mathrm{Equations}12a\text{-}12f\right)\end{array}$
where P_ij is termed the production, Π_ij is termed the pressure-strain correlation tensor, J_ijk is termed the diffusive flux of the Reynolds stresses, and ε_ij is termed the dissipation tensor.

In the general case, equations for all six Reynolds stress terms, and for the energy dissipation rate must expressed. The equation for the isotropic dissipation rate, ε, is similar in structure to the equations for the transport of Reynolds stresses. Full mathematical expressions for the Reynolds stress and isotropic dissipation rate equations shall be expressed in the following section for the specific case of a two-dimensional turbulent boundary layer.

In summary, the equations which are solved to determine turbulent boundary layer parameters include:

• • the continuity equation, Equation 5,
  • the three equations of motion, Equations 9a, 9b, 9c,
  • six equations for Reynolds shear and normal stresses, Equations 12a-12f, and
  • the equation for isotropic dissipation rate, ε (Equation 16 below).

The methodology for the solution of turbulent flow parameters involves a finite difference approximation of the system of equations of motion and continuity, with accompanying boundary conditions.

Turbulent Boundary Layer Equations: Complete mathematical formulations are provided for the specific case of a turbulent boundary layer with a steady, two-dimensional mean flow and a constant freestream velocity, U_∞. Two-dimensional turbulent boundary layer equations, as termed in the literature, are derived from the general continuity equation (Equation 5) and equations of motion (Equations 9a-9c), given the assumptions that:

• • The mean transverse velocity, W, is zero.
  • Gravitational forces can be neglected.
  • The pressure gradient is approximately zero in the y direction.
  • The mean velocity in the streamwise direction, U, is much greater than the mean velocity in the normal direction, V.
  • The rate of change of parameters in the x direction is much smaller than the rate of change of parameters in the y direction.
    The above assumptions permit simplification of the set of equations required to solve for turbulent boundary layer parameters to include the revised continuity equation, Equation (13): $\begin{array}{cc}\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0& \left(\mathrm{Equation}13\right)\end{array}$
    and the equation of motion for the U velocity component (Equation 14): $\begin{array}{cc}U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}=-\frac{1}{\rho }\frac{\partial 𝒫}{\partial x}+v\frac{{\partial }^{2}U}{\partial y^2}-\frac{\partial \stackrel{\_}{u^{\prime }{v}^{\prime }}}{\partial y}& \left(\mathrm{Equation}14\right)\end{array}$

Where transport equations for the six Reynolds stress components are required in the general case, the Reynolds shear stress components − v′w′ and − u′w′ are considered to be very small, so that only equations for u′², v′², w′², and − u′v′ are formulated in the format of Equation 12 (repeated below as Equations 15a-15d). $\begin{array}{cc}{U}_k\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}=P_{\mathrm{ij}}-{\Pi }_{\mathrm{ij}}-\frac{\partial J_{\mathrm{ijk}}}{\partial x_k}-2{\varepsilon }_{\mathrm{ij}}& \left(\mathrm{Equations}15a\text{-}15d\right)\end{array}$
where P_ij is the production term, Π_ij is the pressure-strain correlation tensor, J_ijk is the diffusive flux of the Reynolds stresses, and ε_ij is the dissipation tensor. A fifth equation for ε is: $\begin{array}{cc}U\frac{\partial \varepsilon }{\partial x}+V\frac{\partial \varepsilon }{\partial y}=C_{s1}{f}_1\frac{\varepsilon }{k}{P}_{\Sigma }-C_{s2}{f}_2\frac{\varepsilon }{k}\left[\varepsilon -v\frac{{\partial }^{2}k}{\partial y^2}\right]+\frac{\partial }{\partial y}\left(2{\varepsilon }_t\frac{\partial \varepsilon }{\partial y}\right)+v\frac{{\partial }^2\varepsilon }{\partial y^2}& \left(\mathrm{Equation}16\right)\end{array}$
where the expression for viscous diffusion may alternatively be approximated as: $\begin{array}{cc}[[\hspace{1em}\hspace{1em}v\frac{{\partial }^{2}k}{\partial y^2}=2v\left(\frac{\partial k^{1/2}}{\partial y}\right)& \left(\mathrm{Equation}17\right)\hspace{1em}]\hspace{1em}]\\ v\frac{{\partial }^{2}k}{\partial y^2}\approx 2v\left(\frac{\partial k^{1/2}}{\partial y}\right)& \left(\mathrm{Equation}17\right)\end{array}$
if required for numerical stability in solutions of viscoelastic, non-oscillating surfaces with limited grid points in the near-wall region.

In Equations (15a-15d), the term P_ij may be expressed as: $\begin{array}{cc}{P}_{\mathrm{ij}}=-\stackrel{\_}{u_i^{\prime }{u}_k^{\prime }}\frac{\partial U_j}{\partial x_k}-\stackrel{\_}{u_j^{\prime }{u}_k^{\prime }}\frac{\partial U_i}{\partial x_k}& \left(\mathrm{Equations}18a\text{-}18d\right)\end{array}$
In Equation (16), the term P_Σ may be expressed as: $\begin{array}{cc}{P}_{\Sigma }=\frac{1}{2}{P}_{\mathrm{ii}}& \left(\mathrm{Equation}19\right)\end{array}$

The pressure-strain correlation tensor, Π_ij, which redistributes energy between different components of Reynolds stresses, may be expressed as: $\begin{array}{cc}{\Pi }_{\mathrm{ij}}=C_1\left(\frac{\varepsilon }{k}\right)\left(\stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}-{\delta }_{\mathrm{ij}}\frac{2}{3}k\right)+C_2\left(P_{\mathrm{ij}}-{\delta }_{\mathrm{ij}}\frac{2}{3}{P}_{\Sigma }\right)+{\pi }_{\mathrm{ij},1}^{\prime }+{\pi }_{\mathrm{ij},2}^{\prime }+{\pi }_{\mathrm{ij},3}^{\prime }& \left(\mathrm{Equations}20a\text{-}20d\right)\end{array}$
where the π′_ij,1 terms represent near-wall redistribution of turbulent energy from the streamwise component to the normal and transverse components, the π′_ij,2 terms represent near-wall variation of the Reynolds stress tensor component production, and the π′_ij,3 terms represent near-wall redistribution of turbulent energy proportional to local vorticity: $\begin{array}{cc}\begin{array}{cc}{\pi }_{\mathrm{ij},1}^{\prime }=-C_1^{\prime }\frac{\varepsilon }{k}\left(\stackrel{\_}{v^{\mathrm{\prime 2}}}{\delta }_{ij}-\frac{3}{2}\left(\stackrel{\_}{v^{\prime }{u}_i^{\prime }}{\delta }_{\mathrm{jl}}+\stackrel{\_}{v^{\prime }{u}_j^{\prime }}{\delta }_{\mathrm{il}}\right)\right)f\left(\frac{l}{y}\right)& \end{array}& \left(\mathrm{Equations}21a\text{-}21d\right)\\ {\pi }_{\mathrm{ij},2}^{\prime }=-C_2^{\prime }\left(P_{\mathrm{ij}}-\frac{2}{3}{\delta }_{ij}{P}_{\Sigma }\right)f\left(\frac{l}{y}\right)& \left(\mathrm{Equations}22a\text{-}22d\right)\\ [[\hspace{1em}\hspace{1em}{\pi }_{\mathrm{ij},3}^{\prime }=-C_3^{\prime }\left(P_{\mathrm{ij}}-D_{\mathrm{ij}}\right)f\left(\frac{l}{y}\right)& \left(\mathrm{Equations}23a\text{-}23d\right)\hspace{1em}]\hspace{1em}]\\ {\pi }_{\mathrm{ij},3}^{\prime }=C_3^{\prime }\left(P_{\mathrm{ij}}-D_{\mathrm{ij}}\right)f\left(\frac{l}{y}\right)& \left(\mathrm{Equations}23a\text{-}23d\right)\end{array}$
Here $f\left(\frac{l}{y}\right)$
is a unique damping function for the near-wall region: $\begin{array}{cc}f\left(\frac{l}{y}\right)=\frac{R_t}{R_k}\left(1+\sqrt{1+\frac{100}{R_t}}\right)& \left(\mathrm{Equation}24\right)\end{array}$
where: $\begin{array}{cc}{R}_k=\frac{k^{1/2}y}{v}& \left(\mathrm{Equation}25\right)\end{array}$
and $\begin{array}{cc}{R}_t=\frac{k^2}{v\varepsilon }& \left(\mathrm{Equation}26\right)\end{array}$
Here, D_ij is a dissipation tensor: $\begin{array}{cc}\begin{array}{c}{D}_{\mathrm{ij}}=-\stackrel{\_}{u_i{u}_j}\frac{\partial U_i}{\partial x_j}-\stackrel{\_}{u_j{u}_i}\frac{\partial U_i}{\partial x_i}\\ \frac{\partial J_{\mathrm{ijk}}}{\partial x_k}\end{array}& \left(\mathrm{Equations}27a\text{-}27d\right)\end{array}$
is the gradient of turbulent and viscous diffusive flux of the Reynolds stresses in the boundary layer, where only one component remains in the boundary-layer representation: $\begin{array}{cc}-\frac{\partial J_{\mathrm{ijk}}}{\partial x_2}\pm \underset{\mathrm{turbulent}\mathrm{diffusion}}{\underbrace{\frac{\partial }{\partial y}\left({\left(AC\right)}_t\frac{k}{\varepsilon }\left[\stackrel{\_}{v^{\mathrm{\prime 2}}}\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial y}\right]\right)}}+\underset{\mathrm{viscous}\mathrm{diffusion}}{v\underbrace{\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial y^2}}}& \left(\mathrm{Equations}28a\text{-}28d\right)\end{array}$
where A is 6 in the equation for v′², 2 for u′² and w′², and 4 for − u′v′ (i.e., the effective gradients of turbulent diffusion are different for different components of Reynolds stress), and where the coefficient of turbulent diffusion, ε_t, is: $\begin{array}{cc}{\varepsilon }_t=C_t\frac{k}{\varepsilon }\stackrel{\_}{v^{\mathrm{\prime 2}}}& \left(\mathrm{Equation}29\right)\end{array}$
except for Equation (16), where:
C_i=C_ε  (Equation 30)

The dissipation tensor, ε_ij, is written as: $\begin{array}{cc}{\varepsilon }_{ij}=v\stackrel{\_}{\frac{\partial \stackrel{\_}{u_i^{\prime }}}{\partial x_k}\frac{\partial \stackrel{\_}{u_j^{\prime }}}{\partial x_k}}\approx f_s\frac{\stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{2k}\varepsilon +\left(1-f_s\right)\frac{1}{3}{\delta }_{ij}\varepsilon & \left(\mathrm{Equations}31a\text{-}31d\right)\end{array}$
where ƒ_s characterizes flow in the near-wall region: $\begin{array}{cc}{f}_s=\frac{1}{1+0.06R_t}\mathrm{and}& \left(\mathrm{Equation}32\right)\\ R_t=\frac{k^2}{v\varepsilon }& \left(\mathrm{Equation}33\right)\end{array}$

Equation (16) includes two functions, ƒ₁ and ƒ₂, which also introduce corrections for near-wall flows:
ƒ₁=1+0.8e^−Rt  (Equation 34)
ƒ₂=1−0.2e^−Rt2  (Equation 35)
Values of constants for flow over a flat plate are as shown in Table 1:

C1	C2	Cε1	Ca2	Ct	Cε	C1′	C2′	C3′
1.34	0.8	1.45	1.9	0.12	0.15	0.36	0.45	0.036

For the case of a two-dimensional boundary layer, turbulent boundary layer parameters at different x and y locations are determined through solution of the continuity equation (Equation 13), the equation of motion in the x-direction (Equation 14), the transport equations for the Reynolds stresses u′², v′², w′², and − u′v′ (Equations 15a-15d), and the equation for the isotropic energy dissipation rate (Equation 16), given appropriate boundary conditions. The problem is solved numerically using a finite difference approximation. All equations are reduced to a standard type of parabolic equation in terms of a given function, and solution is obtained at designated grid points in an (x,y) coordinate system.

Boundary Conditions: Boundary conditions are values of parameters at the limits of the boundary layer, i.e., at the surface and the freestream. The freestream velocity is defined as U_∞. Boundary conditions at the surface are specified for Reynolds normal and shear stresses (kinematic boundary conditions), as well as for the isotropic dissipation rate (dynamic boundary condition). For an arbitrary geometry, the x and y coordinates of the surface must be specified. If the surface is a flat plate, the boundary will be along the line y=0.

Since oscillation amplitudes at the surface are small, linearized kinematic boundary conditions, where mean velocities at the surface are assumed to be zero, are appropriate. Boundary conditions for fluctuating velocity components at the surface of a flat plate are expressed as: $\begin{array}{cc}{u^{\prime }}_{y=0}=\frac{\partial {\xi }_1}{\partial t}\cos \Theta -{\xi }_2\frac{u_{*}^2}{v}& \left(\mathrm{Equation}36\right)\\ {v^{\prime }}_{y=0}=\frac{\partial {\xi }_2}{\partial t}& \left(\mathrm{Equation}37\right)\\ {w^{\prime }}_{y=0}=\frac{\partial {\xi }_1}{\partial t}\sin \Theta & \left(\mathrm{Equation}38\right)\end{array}$
where ξ₁ and ξ₂ are the longitudinal and vertical surface displacement components, respectively, u* is the friction velocity (as previously defined), and Θ is the angle of the longitudinal axis relative to the mean flow in the x₁-x₃ plane. With linearized boundary conditions, mean velocities at the wall are assumed to be zero. Surface displacements are approximated by the first mode of a Fourier series: $\begin{array}{cc}{\xi }_i=\sum _{j=1}^{\infty }{a}_{\mathrm{ij}}{e}^{i{\alpha }_j\left(x-\mathrm{Ct}\right)}\approx a_{\mathrm{il}}{e}^{i{\alpha }_e\left(x-\mathrm{Ct}\right)}& \left(\mathrm{Equation}39\right)\end{array}$
Here, α₃ is the wavenumber corresponding to the maximum turbulent energy in the boundary layer, and is given by: $\begin{array}{cc}{\alpha }_e=\frac{{\omega }_e}{C}& \left(\mathrm{Equation}40\right)\end{array}$
where the energy-carrying frequency, ω_e, is assumed to be: $\begin{array}{cc}{\omega }_e=\frac{U_{\infty }}{\delta }& \left(\mathrm{Equation}41\right)\end{array}$
and the phase speed corresponding to energy-carrying disturbances in the boundary layer is assumed to be:
C≈0.8U_∞  (Equation 42)
Since there is a range of frequencies which carry energy, as reported within the scientific literature, it is advantageous to also perform calculations for the case where: $\begin{array}{cc}{\omega }_e=2\frac{U_{\infty }}{\delta }& \left(\mathrm{Equation}43\right)\end{array}$
In the absence of resonance, it is appropriate to time-average components of the Reynolds stress at the wall: $\begin{array}{cc}{\stackrel{\_}{u^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2\left({{\xi }_1}^2+\frac{2u_{*}^2}{{\omega }_{e}v}{\xi }_1{\xi }_2\sin \left({\phi }_2-{\phi }_1\right)+\frac{1}{{\omega }_e^2}{\left(\frac{u_{*}^2}{v}\right)}^2{{\xi }_2}^2\right)& \left(\mathrm{Equation}44\right)\\ {\stackrel{\_}{v^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2{{\xi }_2}^2& \left(\mathrm{Equation}45\right)\\ {\stackrel{\_}{w^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2{{\xi }_1}^2{\tan }^2\Theta & \left(\mathrm{Equation}46\right)\\ {-\stackrel{\_}{u^{\prime }{v}^{\prime }}}_{y=0}=\frac{1}{2}{\omega }_e^2{\xi }_2{\xi }_1\cos \left({\phi }_2-{\phi }_1\right)& \left(\mathrm{Equation}47\right)\end{array}$
where |ξ_i| is the rms amplitude of the displacement. For a passive isotropic viscoelastic coating excited by a forced load, the response takes the form of a traveling wave, so that the phase shift between normal and longitudinal
where ε_ij^s is the strain tensor: $\begin{array}{cc}{\varepsilon }_{\mathrm{ij}}^s=\frac{1}{2}\left(\frac{\partial {\xi }_i}{\partial x_j}+\frac{\partial {\xi }_j}{\partial x_i}\right)& \left(\mathrm{Equations}53a\text{-}53d\right)\end{array}$
and:
ε^s=ε_ii^s  (Equation 54)
λ(ω) is the frequency-dependent Lame constant, which is defined in terms of the bulk modulus, K(ω), which can be reasonably approximated as the static bulk modulus, K₀, and the complex shear modulus, μ(ω): $\begin{array}{cc}\lambda \left(\omega \right)=K_0-\frac{2}{3}\mu \left(\omega \right)& \left(\mathrm{Equation}55\right)\end{array}$

Displacements, ξ_i, are approximated as periodic, in the form of Equation (39), and can be expressed as a function of potentials of longitudinal and transverse (shear) waves: $\begin{array}{cc}\stackrel{\_}{\xi }=\nabla \varphi +\nabla \overrightarrow{\Psi }\overrightarrow{\Psi }=\left\{0,0,\Psi \right\}& \left(\mathrm{Equation}56\right)\end{array}$
where ∇_φ is the gradient of φ and ∇×{right arrow over (Ψ)} is the curl of the vector {right arrow over (Ψ)}. Equation (54) can be rewritten as two decoupled equations for the two wave potentials: $\begin{array}{cc}\left[2\mu \left(\omega \right)+\lambda \left(\omega \right)\right]{\nabla }^2\phi ={\rho }_s\frac{{\partial }^2\phi }{\partial t^2}& \left(\mathrm{Equation}57\right)\\ \mu \left(\omega \right){\nabla }^2\Psi ={\rho }_s\frac{{\partial }^2\Psi }{\partial t^2}& \left(\mathrm{Equation}58\right)\end{array}$

Equations (57) and (58) can be solved for the potentials, φ and Ψ, and hence for displacements, velocities, and stresses through the thickness of the coating, if boundary conditions are specified. The coating is fixed at its base, so that the longitudinal and normal displacements are zero, and the shear stress and pressure load on the surface is known. Pressure and shear pulsations on the coatings are approximated as periodic functions, with a form similar to that of the displacements in Equation (39), but with the following magnitudes, respectively:




where K_p is the Kraichnan parameter, whose value is approximated as 2.5.

If shear pulsations are included, a phase shift between shear and pressure pulsations must also be introduced.

If calculations are performed for a unit load, then surface displacements under actual load will be: $\begin{array}{cc}{\xi }_i{❘}_{p_{\mathrm{rms}=\mathrm{actual}}}=\frac{\rho K_p{u}_{*}^2{\xi }_i{❘}_{p_{\mathrm{rms}}=1}}{H}H{\tilde{u}}_{*}^2=C_{\mathrm{ki}}H{\tilde{u}}_{*}^2& \left(\mathrm{Equation}61\right)\end{array}$
Kinematic boundary conditions in Equations (44) to (47) for the turbulent boundary layer problem are rewritten in terms of output from the materials problem: $\begin{array}{cc}\frac{\stackrel{\_}{u^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{\left(C_{\mathrm{k1}}+{\tilde{u}}_{*}{\mathrm{Re}}_{*}{C}_{\mathrm{k2}}\right)}^2{\tilde{u}}_{*}^4\left({\mathrm{Re}}_{*}=\frac{u_{*}\delta }{v}\right)& \left(\mathrm{Equation}62\right)\\ \frac{\stackrel{\_}{v^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{C}_{\mathrm{k2}}^2{\tilde{u}}_{*}^4& \left(\mathrm{Equation}63\right)\\ \frac{\stackrel{\_}{w^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{C}_{\mathrm{k1}}^2{\tilde{u}}_{*}^4{\tan }^2\Theta & \left(\mathrm{Equation}64\right)\\ -\frac{\stackrel{\_}{u^{\prime }{v}^{\prime }}}{U_{\infty }^2}=0& \left(\mathrm{Equation}65\right)\end{array}$

Dynamic boundary conditions are rewritten in the form: $\begin{array}{cc}{{-\frac{{\stackrel{\_}{p^{\prime }{v}^{\prime }}}_{y=0}}{{\mathrm{\rho U}}_{\infty }^3}\approx {\tilde{\varepsilon }}_t\frac{\partial \tilde{k}}{\partial \tilde{y}}}_{}}_{y=0}& \left(\mathrm{Equation}66\right)\end{array}$
where:
C_k3=C_k2K_pγ(ω)
and where γ(ω) is a dissipative function of the coating material.

The flux of turbulent fluctuating energy through the surface can be solved for directly, as:
− p′v′|_y=0=−¼|p′|ω(−tξ₂+tξ*₂)  (Equation 67)
but the nondimensionalized flux can also be approximated as a diffusive flux term, using the gradient diffusion approach: $\begin{array}{cc}{{-\frac{{\stackrel{\_}{p^{\prime }{v}^{\prime }}}_{y=0}}{{\mathrm{\rho U}}_{\infty }^3}\approx {\tilde{\varepsilon }}_t\frac{\partial \tilde{k}}{\partial \tilde{y}}}_{}}_{y=0}& \left(\mathrm{Equation}68\right)\end{array}$
where: $\begin{array}{cc}\tilde{k}=\frac{k}{U_{\infty }^2}& \left(\mathrm{Equation}69\right)\\ \tilde{y}=\frac{y}{\delta }& \left(\mathrm{Equation}70\right)\\ {\tilde{\varepsilon }}_t=\frac{{\varepsilon }_t}{U_{\infty }\delta }& \left(\mathrm{Equation}71\right)\end{array}$

Equation 68 provides a basis for determining the value of the kinematic coefficient of turbulence diffusion, ε_q, on an absorbing surface. Substituting Equation 68 into Equation
based on the criterion that the speed of shear waves in the material is approximately the same as the phase speed of the energy-carrying disturbances, C. This phase speed, C, is assumed to be 0.8 of the value of the freestream velocity, U_∞ (Equation 42). If the convective velocity exceeds the shear wave velocity, an instability occurs, and large waves appear on the surface of the material, leading to an increase of drag for the coating.

An initial choice for thickness, H, for isotropic viscous materials, where $\frac{{\mu }_s}{{\mu }_0}>1$
is: $\begin{array}{cc}H=\frac{3C}{{\omega }_e}& \left(\mathrm{Equation}79\right)\end{array}$
and for isotropic, low viscosity materials, where $\frac{{\mu }_s}{{\mu }_0}<1$
is: $\begin{array}{cc}H=\frac{5C}{{\omega }_e}& \left(\mathrm{Equation}80\right)\end{array}$

The optimal desired thickness for a coating may be greater than practical for a given application. While isotropic coatings thinner than recommended in Equations 79-80 can still be effective, anisotropic coatings that are stiffer in the normal dimension relative to the transverse and longitudinal dimensions can provide equivalent performance with significant reduction in thickness.

Given specified values of H, μ₀, and ρ_s, the VE problem, as expressed in Equations (57) and (58), is solved numerically for a matrix of values of τ_s and μ_s (i.e., for different values of the complex shear modulus) and for a range of wavenumbers. The wavenumber corresponding to the maximum turbulent energy in the boundary layer is: $\begin{array}{cc}{\alpha }_e=\frac{{\omega }_e}{C}& \left(\mathrm{Equation}81\right)\end{array}$
where the frequency, ω_e, for maximum energy-carrying disturbances is estimated by Equation (41). Calculations yield surface displacement amplitudes and the flux of turbulent fluctuating energy into the coating. The best combination of properties for a SRT material occurs where the surface displacement under actual load (Equation (61)) is less than the viscous sublayer thickness, and where the energy flux into the coating (Equation (57)) is at a maximum. Furthermore, it is desirable to maintain this criterion for a range of frequencies from approximately one decade below to one decade above the energy-carrying frequency, ω_e.

Once a set of optimal values of τ_s and μ_s are determined for a set of specified values of H, μ₀, and ρ_s, the calculations are iterated using slightly different values of thickness, H, and static modulus, μ₀. From these calculations are chosen the optimal set of parameters for an SRT material (H, μ₀, ρ_s, τ_s and μ_s), given specified flow conditions and configuration.

The complex shear moduli of real polymeric materials, such as polyurethanes and silicones, which are candidates for viscoelastic coatings cannot be adequately described by the SRT representation. More complex MRT or HN representations of the shear modulus require multiple constants, and are less suitable for numerical parametric evaluation. Therefore, results for SRT materials are used to select candidate materials, such as described by the HN formulation (Equation (79)) which can be more readily fabricated in practice. As a guideline, it is desired to match the complex shear modulus curves of the target SRT material and the HN material (value and slope) over frequencies ranging from decade below to one decade above ω_e, with the most important matching being in the immediate vicinity of ω_e.

To design a multi-layer isotropic coating, properties of a complex shear modulus, density, and thickness are specified for individual layers, and non-slip boundary conditions between layers are imposed. The properties of the upper layer are specified according to the methodology for a single layer, and the lower layers will have progressively lower static shear moduli, as optimized for lower freestream velocities. Thus, well-designed multi-layer coatings can reduce drag over a range of freestream velocities.

In the design of an anisotropic coating, the complex shear modulus has different values in the normal direction relative to the longitudinal and transverse directions (hereinafter termed transversely isotropic). If the viscoelastic material follows a single-relaxation time model, then the static shear modulus, μ₀, the dynamic shear modulus, μ_s, and the relaxation time, τ_s, will differ with direction, as expressed in Equations (82) and (83). The static shear modulus in the normal direction, μ₀₁, will be greater than than in the longitudinal-transverse plane, μ₀₂. The complex shear modulus in the normal direction is expressed as: $\begin{array}{cc}{\mu }_1\left(\omega \right)={\mu }_{01}+{\mu }_{\mathrm{s1}}\left[\frac{{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}+i\frac{{\mathrm{\omega \tau }}_{\mathrm{s1}}}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}\right]& \left(\mathrm{Equation}82\right)\end{array}$
while the shear modulus in the streamwise and transverse directions is expressed as: $\begin{array}{cc}{\mu }_2\left(\omega \right)={\mu }_{02}+{\mu }_{\mathrm{s2}}\left[\frac{{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}+i\frac{{\mathrm{\omega \tau }}_{\mathrm{s2}}}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}\right]& \left(\mathrm{Equation}83\right)\end{array}$
For a viscoelastic, transversely isotropic material, surface oscillation amplitudes can be reduced relative to an isotropic material, while the level of energy flux into the material is increased. Thus, well-designed anisotropic coatings will be significantly thinner than isotropic coatings associated with the same level of drag reduction.
Methodology to Choose Structure of a Drag-Reducing Viscoelastic Coating

A further aspect of coating design is the choice of internal structure within the viscoelastic material. In practical applications of viscoelastic coatings, the coating will be finite in length, with leading, trailing, and side edges. The influence of the finite edges affects coating performance. Well posed edges can order and stabilize transverse and longitudinal vortical structures in the near-wall region of the flow and thereby delay the deformation of these vortical structures and enhance the stability of the flow. However, unstructured edges can accentuate the amplitude of oscillations of the viscoelastic material in this region. Local instabilities can degrade the performance of the coating, so that, even with a well-designed material, the influence of the edges can lead to a drag increase. Hence, the coating is structured in the vicinity of finite edges. The thickness of the coating is decreased to minimize such oscillations, using techniques such as a rigid wedge underneath the coating, or other localized structure near an edge (FIG. 3). For large bodies, a continuous coating may be impractical or difficult to fabricate. An alternative design is a piecewise continuous coating, composed of finite segments of coating, where both the longitudinal and transverse edges of the coating system are organized to stabilize flow structures and to minimize adverse effects at the edges of each segment.

In addition to well-posed edges, viscoelastic coatings may be combined with surface structure to enhance the stabilization of longitudinal vortices along the length of the coating, and hence to increase the level of drag reduction through multiple physical mechanisms. Structure can include the placement of riblets on top of the viscoelastic coating, or the creation of so-called “inverse” riblets. In the latter case, a viscoelastic coating may be molded over ribs or ridges of rigid material, so that longitudinal riblet structures form when fluid flows over the viscoelastic surface.

The dimensions (scales) of the segments and the dimensions of the structures within the coating are selected as multiples of the transverse and longitudinal scales in the near wall turbulent flow. These scales vary with body speed, position along the body and when non-Newtonian additives, such as dilute aqueous solutions of high-molecular weight polymers, are present.

Claims

1. A method for estimating the reduction in friction drag of a body with a viscoelastic coating as compared to the friction drag of a rigid surface of identical size and shape as said body, during turbulent flow at a specified freestream velocity, said method comprising the following steps, performed in the order indicated:

(I) using boundary conditions for a rigid surface, determining characteristics of a turbulent boundary layer at the specified freestream velocity, said characteristics including the boundary layer thickness, phase speed and frequency corresponding to maximum-energy-carrying disturbances, mean velocity profiles, Reynolds stress distributions, wall shear stress distribution, and friction drag,

(II) selecting a density, complex shear modulus, and thickness of a viscoelastic coating which corresponds to minimum oscillation amplitudes and maximum flux of energy into the viscoelastic coating during excitation by a forcing function substantially identical to the excitation produced by the turbulent boundary layer as determined in step I, and

(III) using oscillation amplitudes and energy flux corresponding to the viscoelastic coating of selected density, complex shear modulus and thickness as determined in step II, determining characteristics of the turbulent boundary layer at the specified freestream velocity, including mean velocity profiles, Reynolds stress distributions, wall shear stress distribution, and friction drag, and

(IV) determining a percentage reduction in friction drag as the ratio of i) the friction drag as determined in step I minus the friction drag as determined in step III, divided by ii) the friction drag as determined in step I.

2. The method of claim 1, wherein step I thereof includes the following substeps, not necessarily performed in the order indicated:

(a) for a rigid surface of specified geometry, and for a given freestream velocity, U_∞, solving the equation of continuity: $\frac{\partial U_i}{\partial x_i}=0$

 and the equations of motion for an incompressible fluid and steady flow with constant kinematic viscosity, ν, constant density, ρ, negligible body forces, and gradients of pressure, P: $U_j\frac{\partial U_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial 𝒫}{\partial x_i}+g_i+v\frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}$

where U_i refers to velocity components in the x, y, and z directions (x₁, x₂, and x₃ in indicial notation) and where fluid velocities at the surface of the body are zero,

(b) from the velocity field, U_i, determined from the solution of the general continuity and motion equations in step 2(a), determining the boundary layer thickness as a function of body geometry, where the edge of the boundary layer is defined as that location where the ratio of mean velocity to the freestream velocity is a constant, β, between 0.95 and 1.0: ${\frac{U}{U_{\infty }}}_{\frac{y}{\delta }=1}=\beta$

(c) from the solution of the general continuity and motion equations in step 2(a), determining the shear stress, τ_w, along the body, the local coefficient of friction drag, and the integral drag, and

(d) estimating the frequency, ω_e, and the phase speed, C, corresponding to maximum-energy-carrying disturbances in the boundary layer.

3. The method of claim 2, wherein in substep 2(b), the boundary layer thickness is approximated as that location where the ratio of the mean velocity to the freestream velocity has a constant value of β=9975. velocity has a constant value of β=0.9975.

4. The method of claim 2, wherein in substep 2(d), the maximum energy-carrying frequency for disturbances in the boundary layer is estimated as: ${\omega }_e=\frac{U_{\infty }}{\delta }$

5. The method of claim (2), for the case of two-dimensional flow over a flat plate, wherein the equations of continuity and the equations of motion are reduced to the following system of equations, where a Reynolds transport approach is adopted for turbulence closure, and where P is the mean pressure, ν is the kinematic viscosity, ρ is the density, U is the mean longitudinal velocity component, V is the mean normal velocity component, and u′, v′, and w′ are components of fluctuating velocity in the streamwise, normal, and transverse directions, respectively, and ε is the isotropic dissipation rate: $\begin{array}{cc}\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0& \\ U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}=-\frac{1}{\rho }\frac{\partial 𝒫}{\partial x_i}+v\frac{{\partial }^{2}U}{\partial y^2}-\frac{\partial \stackrel{\_}{u^{\prime }{v}^{\prime }}}{\partial y}& \\ U_k\frac{\partial \stackrel{\_}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}=P_{\mathrm{ij}}-{\Pi }_{\mathrm{ij}}-\frac{\partial J_{\mathrm{ijk}}}{\partial x_k}-2{\varepsilon }_{\mathrm{ij}}& \\ U\frac{\partial \varepsilon }{\partial x}+V\frac{\partial \varepsilon }{\partial y}=C_{\varepsilon 1}{f}_1\frac{\varepsilon }{k}{P}_{\Sigma }-C_{\varepsilon 2}{f}_2\frac{\varepsilon }{k}\left[\varepsilon -v\frac{{\partial }^{2}k}{\partial y^2}\right]+\frac{\partial }{\partial y}\left(2{\varepsilon }_t\frac{\partial \varepsilon }{\partial y}\right)+v\frac{{\partial }^2\varepsilon }{\partial y^2}& \end{array}$

6. The method of claim 5 wherein Reynolds stress boundary conditions are derived from values of surface amplitudes, ξ₁ and ξ₂, calculated from solution of the momentum equation for a viscoelastic material: $\begin{array}{c}{\stackrel{\_}{u^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2\left({{\xi }_1}^2+\frac{2u_{*}^2}{{\omega }_{e}v}{\xi }_1{\xi }_2\sin \left({\phi }_2-{\phi }_1\right)+\frac{1}{{\omega }_e^2}{\left(\frac{u_{*}^2}{v}\right)}^2{{\xi }_2}^2\right)\\ {\stackrel{\_}{v^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2{{\xi }_2}^2\\ {\stackrel{\_}{w^{\mathrm{\prime 2}}}}_{y=0}=\frac{1}{2}{\omega }_e^2{{\xi }_1}^2{\tan }^2\Theta \\ {-\stackrel{\_}{u^{\prime }{v}^{\prime }}}_{y=0}=\frac{1}{2}{\omega }_e^2{\xi }_2{\xi }_1\cos \left({\phi }_2-{\phi }_1\right)\end{array}$

7. The method of claim 6 wherein the energy boundary condition is expressed in terms of the solution of the momentum equation for a viscoelastic material where:

(a) the energy flux into the coating is determined from the pressure velocity correlation, which, given a periodic forcing function can be expressed in terms of the amplitude of the normal oscillation amplitude (and its complex conjugate, as designated by *): ${-\stackrel{\_}{p^{\prime }{v}^{\prime }}}_{y=0}=-\frac{1}{4}p^{\prime }\omega \left(-{\mathrm{i\xi }}_2+{\mathrm{i\xi }}_2^{*}\right)$

(b) the amplitude of surface oscillations, ξ₂, given surface loading corresponding to the turbulent boundary layer, is expressed as: ${\xi }_2{❘}_{p_{\mathrm{rms}=\mathrm{actual}}}=\frac{\rho K_p{u}_{*}^2{\xi }_2{❘}_{p_{\mathrm{rms}}=1}}{H}H{\tilde{u}}_{*}^2=C_{\mathrm{k2}}H{\tilde{u}}_{*}^2$

 where:

H is the thickness of the coating, ${\tilde{u}}_{*}=\frac{u_{*}}{U_{\infty }},\mathrm{and}$

u* is the friction velocity, defined as μ₀ is the friction velocity, defined as $u_{*}=\sqrt{\frac{{\tau }_w}{\rho }},$

 where:

ρ is the density of the fluid, and

τ_w is the shear stress at the wall, expressed as ${{\tau }_w=\rho v\frac{\partial U}{\partial y}}_{y=0}$

 so that the pressure velocity correlation can be expressed as: $-\frac{{-\stackrel{\_}{p^{\prime }{v}^{\prime }}}_{y=0}}{{\mathrm{\rho U}}_{\infty }^3}=-\frac{1}{4}\frac{H}{\delta }{C}_{\mathrm{k3}}{\stackrel{\_}{u}}_{*}^4$

where:

C_k3=C_k2K_pγ(ω)

and where K_p is the Kraichnan parameter, with a value of 2.5, and γ(ω) is a dissipative function of the coating material, which reflects the phase shift between the pressure fluctuations and the vertical displacement of the coating,

(c) the flux of energy into the coating can be approximated by an effective turbulent diffusion term, ${\varepsilon }_t\frac{\partial k}{\partial y},$

 where: ${\varepsilon }_t=C_t\frac{k}{\varepsilon }\stackrel{\_}{v^{\mathrm{\prime 2}}},$

 where C_t=0.12

 so that, in terms of nondimensionalized quantities: ${-\frac{{-\stackrel{\_}{p^{\prime }{v}^{\prime }}}_{y=0}}{{\mathrm{\rho U}}_{\infty }^3}\approx {\tilde{\varepsilon }}_t\frac{\partial \tilde{k}}{\partial \tilde{y}}}_{y=0}$

 where: $\begin{array}{c}{\tilde{\varepsilon }}_t=\frac{{\varepsilon }_t}{U_{\infty }\delta }\\ \tilde{k}=\frac{k}{U_{\infty }^2}\end{array}$

e####

and: $\tilde{y}=\frac{y}{\delta }$

(d) given the expression for nondimensionalized pressure-velocity correlation in step (b) and in step (c), ε_q=ε_t|_y=0 can be written as: ${\tilde{\varepsilon }}_q=\frac{{\varepsilon }_t{❘}_{y=0}}{U_{\infty }\delta }=\frac{1}{4}{C}_{\mathrm{k3}}\left(\frac{H}{\delta }\right){\tilde{u}}_{*}\frac{{\tilde{u}}_{*}^3{\Delta }_{\max }}{k_{\max }^{+}-k_q^{+}}$

 where γ_max⁺ is the normal distance from the surface to the maximum of turbulence energy, nondimensionalized as follows: $\begin{array}{c}{y}_{\max }^{+}=\frac{y_{\max }}{v}{u}_{*}\\ {\Delta }_{\max }=\frac{y^{+}}{{\mathrm{Re}}_{*}}\end{array}$

where: ${\mathrm{Re}}_{*}=\frac{u_{*}\delta }{v}$

 is the Reynolds number based on friction velocity, u*, and boundary layer thickness, δ

u* is the friction velocity, defined as

μ₀, and boundary layer thickness, δ

μ₀ is the friction velocity, defined above as $u_{*}=\sqrt{\frac{{\tau }_w}{\rho }},$

where:

ρ is the density of the fluid, and

τ_w is the shear stress at the wall, expressed as ${{\tau }_w=\rho v\frac{\partial U}{\partial y}}_{y=0}$

and $k_{\max }^{+}=\frac{k_{\max }}{U_{\infty }^2}$

 is the maximum of turbulence kinetic energy, non-dimensionalized by U_∞² $k_q^{+}=\frac{k_q}{U_{\infty }^2}$

 is the kinetic energy of the oscillating surface, non-dimensionalized by U_∞²

(e) the boundary condition for isotropic dissipation rate is expressed dimensionally as: ${{\varepsilon }_{y=0}=\frac{\partial }{\partial y}\left[\left(v+{\varepsilon }_t\right)\frac{\partial k}{\partial y}\right]}_{y=0}$

 or, in nondimensional form: ${{\tilde{\varepsilon }}_{y=0}=\frac{\partial }{\partial y}\left[\left(\frac{1}{\mathrm{Re}}+{\tilde{\varepsilon }}_t\right)\frac{\partial \tilde{k}}{\partial \tilde{y}}\right]}_{y=0}$

where: $\mathrm{Re}=\frac{U_{\infty }\delta }{v}$

 is the Reynolds number based on freestream velocity and boundary layer thickness.

8. The method of claim 7, where the solution of the turbulent boundary layer problem for an isotropic viscoelastic surface is based upon an iterative technique, and where initial values are assumed for u* and for the gradient of turbulent kinetic energy, ${\frac{\partial k}{\partial y}}_{y=0}$

as based upon the solution of the turbulent boundary layer problem for a rigid flat plate, with two-dimensional flow, where the boundary conditions are: $\begin{array}{c}{\frac{\stackrel{\_}{u^{\mathrm{\prime 2}}}}{U_{\infty }^2}}_{y=0}=0\\ {\frac{\stackrel{\_}{v^{\mathrm{\prime 2}}}}{U_{\infty }^2}}_{y=0}=0\\ {\frac{\stackrel{\_}{w^{\mathrm{\prime 2}}}}{U_{\infty }^2}}_{y=0}=0\\ {\frac{\stackrel{\_}{u^{\prime }{v}^{\prime }}}{U_{\infty }^2}}_{y=0}=0\\ {{\varepsilon }_{y=0}=v\frac{{\partial }^{2}k}{\partial y^2}}_{y=0}\end{array}$

where the values of u* and the gradient of turbulent kinetic where the values of μ₀ and the gradient of turbulent kinetic energy, ${\frac{\partial k}{\partial y}}_{y=0},$

obtained from this solution are used to determine boundary conditions for the next iteration, and where this procedure continues until the solution converges.

9. The method for estimating the reduction in friction drag according to claim 8, for the case of two-dimensional flow over a flat plate coated by an isotropic viscoelastic material, wherein:

a.) the value of the shear modulus, μ(ω), is constant in all directions,

b.) the phase difference between normal and longitudinal displacements at the wall is equal to π/2, so that the Reynolds shear stresses at the surface of the viscoelastic coating are equal to zero:

− u′v′=0

c) Reynolds stress boundary conditions are expressed in non-dimensional form as: $\begin{array}{c}\frac{\stackrel{\_}{u^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{\left(C_{\mathrm{k1}}+{\tilde{u}}_{*}{\mathrm{Re}}_{*}{C}_{\mathrm{k2}}\right)}^2{\tilde{u}}_{*}^4\\ \frac{\stackrel{\_}{v^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{C}_{\mathrm{k2}}^2{\tilde{u}}_{*}^4\\ \frac{\stackrel{\_}{w^{\mathrm{\prime 2}}}}{U_{\infty }^2}=\frac{1}{2}{\left(\frac{H}{\delta }\right)}^2{C}_{\mathrm{k1}}^2{\tilde{u}}_{*}^4{\tan }^2\Theta \\ -\frac{\stackrel{\_}{u^{\prime }{v}^{\prime }}}{U_{\infty }^2}=0\end{array}$

where ${\mathrm{Re}}_{*}=\frac{u_{*}\delta }{v}$

 is the Reynolds number based on friction velocity, u*, and boundary layer thickness, δ

u* is the friction velocity, defined as

μ₀, and boundary layer thickness, δ

μ₀ is the friction velocity, defined above as $u_{*}=\sqrt{\frac{{\tau }_w}{\rho }},$

 where:

ρ is the density of the fluid, and

τ_ωis the shear stress at the wall, expressed as $\begin{array}{c}\begin{array}{c}{{\tau }_w=\rho v\frac{\partial U}{\partial y}}_{y=0}\\ {\tilde{u}}_{*}=\frac{u_{*}}{U_{\infty }}\\ C_{\mathrm{ki}}=\frac{\rho K_p{u}_{*}^2{\xi }_i{❘}_{p_{\mathrm{rms}}=1}}{H}\end{array}\\ K_p=2.5\end{array}$

and where C_k1 and C_k2 can be approximated as zero if the amplitude of surface oscillations is less than the thickness of the viscous sublayer.

10. The method according to claim 8, wherein the coating is composed of multiple layers of isotropic viscoelastic materials with different material properties, where:

a.) the shear modulus, μ(ω), is constant within each of the multiple layers;

b.) boundaries between layers are fixed, and

c.) the static shear modulus of the material is progressively lower for each layer from the top layer to the bottom layer of the coating.

11. The method according to claim 8, wherein the coating is composed of an anisotropic material, where the properties of the material in the normal direction (y) differ from those in the transverse (x-z) plane, wherein the shear modulus in the streamwise and transverse directions is given by μ₁(ω), and the shear modulus in the normal direction is given by μ₂(ω): $\begin{array}{c}{\mu }_1\left(\omega \right)={\mu }_{01}+{\mu }_{\mathrm{s1}}\left[\frac{{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}+i\frac{{\mathrm{\omega \tau }}_{\mathrm{s1}}}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s1}}\right)}^2}\right]\\ {\mu }_2\left(\omega \right)={\mu }_{02}+{\mu }_{\mathrm{s2}}\left[\frac{{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}+i\frac{{\mathrm{\omega \tau }}_{\mathrm{s2}}}{1+{\left({\mathrm{\omega \tau }}_{\mathrm{s2}}\right)}^2}\right].\end{array}$

12. The method of claim 1, wherein step II thereof includes the following substeps, not necessarily performed in the order indicated:

(a) choosing a density, ρ_s, for the coating which is within 10% of the density of water,

(b) selecting an initial choice for the static shear modulus, μ₀, of the material to avoid resonance conditions, based on the criterion that the speed of shear waves in the material, $\sqrt{\frac{{\mu }_0}{{\rho }_s}},$

 is approximately the same as the phase speed, C, of energy-carrying disturbances,

(c) selecting an initial choice for coating thickness, H,

(d) solving the momentum equation for a viscoelastic material: ${\rho }_s\frac{{\partial }^2{\xi }_i}{\partial t^2}=\frac{\partial {\sigma }_{\mathrm{ij}}}{\partial x_j}$

 given harmonic loading of unit amplitude, phase velocity, C, corresponding to the load of the turbulent boundary layer, and variable wavenumber and for a coating which is fixed at its base to a rigid substrate, where:

ξ₁ and ξ₂ are the longitudinal and normal displacements through the thickness of the coating, approximated by the first mode of a Fourier series as follows:

ξ_i=a_i1e^iae(x−Ct)

where a_i1 is a coefficient, α_e=ω_e/C is the wavenumber corresponding to maximum energy in the boundary layer, and φ_i is a phase difference,

σ_ij is the stress tensor, written for a Kelvin-Voigt type viscoelastic material as:

σ_ij=λ(ω)ε^sδ_ij+2μ(ω)ε_ij^s

 where ε_ij^s is the strain tensor, expressed as: ${\varepsilon }_{\mathrm{ij}}^s=\frac{1}{2}\left(\frac{\partial {\xi }_i}{\partial x_j}+\frac{\partial {\xi }_j}{\partial x_i}\right)$

 where ε^s=ε_ii^s,

λ(ω) is the frequency-dependent Lame constant, defined in terms of the bulk modulus, K(ω), and the complex shear modulus, λ(ω): $\lambda \left(\omega \right)=K_0-\frac{2}{3}\mu \left(\omega \right)$

μ(ω) is constant in all directions for an isotropic material, but will have different values in different directions for an anisotropic material,

 and where the momentum equation is solved for a series of materials whose shear moduli can be approximated for a single-relaxation type (SRT) of material, characterized mathematically in the form: $\mu \left(\omega \right)={\mu }_0+{\mu }_s\left[\frac{{\left({\mathrm{\omega \tau }}_s\right)}^2}{1+{\left({\mathrm{\omega \tau }}_s\right)}^2}+i\frac{{\mathrm{\omega \tau }}_s}{1+{\left({\mathrm{\omega \tau }}_s\right)}^2}\right]$

 with an initial value of static shear modulus, μ₀=μ|_ω=0, determined in step 3(b), and for different values of dynamic shear modulus, μ_s, and relaxation time, τ_s, where:

μ_s=μ|_ω=∞−μ₀

(e) from the series of SRT materials for which calculations were performed in step 3(d), selecting a material where the calculated oscillation amplitude does not exceed the viscous sublayer thickness under given flow conditions, and where the energy flux into the coating is considered over a range of frequencies approximately one decade above and below the energy-carrying frequency, ω_e,

(f) iterating steps 3(d) and 3(e) with different values of thickness and static modulus, and determine the optimal combination of material properties for a single relaxation time (SRT) type of material,

(g) using the results from step 3(f) as a guideline in the frequency range of interest, solving the momentum equation for viscoelastic materials whose complex shear modulus may be characterized in the form of a Havriliak-Negami, or HN, material, whose shear modulus is expressed in the following form: $\frac{\mu -{\mu }_{\infty }}{{\mu }_0-{\mu }_{\infty }}=\frac{1}{{\left[1+{\left({\mathrm{i\omega \tau }}_{\mathrm{HN}}\right)}^{{\alpha }_{\mathrm{HN}}}\right]}^{{\beta }_{\mathrm{HN}}}}$

where μ_∞=μ|_ω=∞, τ_HN is a relaxation time, and α_HN and β_HN are constants, and

(h) selecting final material properties of a coating material based on conditions of minimal oscillation amplitudes and maximum energy flux, which is equivalent to the correlation between pressure and velocity fluctuations, − p′v′, in the range of frequencies corresponding to maximum-energy-carrying disturbances, as well as based on consideration of material fabricability.

13. The method of claim 1, wherein step III thereof includes the following substeps, performed in the order indicated:

(a) solving the equations of motion and continuity for a body with a viscoelastic coating, given the same location and flow conditions as in step I, using a numerical methodology which accounts for non-zero energy flux and surface oscillation boundary conditions, as well as for redistribution of energy between fluctuating components in the near-wall region, and

(b) comparing the friction drag calculated for a viscoelastic plate with that calculated for a rigid plate.

14. The method of claim 1, wherein said viscoelastic coating is configured with multiple surface or internal structures that approximate the natural longitudinal and transverse wavelengths in the nearwall flow.

15. The method of claim (1), wherein said viscoelastic coating is structured with a rigid underlying wedge or contour shape to minimize coating thickness and oscillations near the intersection of the coating with a rigid surface.

17. The method according to claim 16, wherein the viscoelastic coating is composed of multiple layers of isotropic, viscoelastic materials with different material properties, where:

a) the complex shear modulus, μ(ω), is constant within each of the multiple layers,

b) boundaries between layers are fixed, and

c) the static shear modulus μ₀ of the material is progressively lower for each layer from the top layer to the bottom layer of the viscoelastic coating.

18. The method according to claim 16, wherein the coating is composed of an anisotropic material, and the properties of the anisotropic material in the normal direction (y) differ from those in the transverse (x-z) plane.

19. The method of claim 16, wherein said viscoelastic coating is structured with a more rigid underlying wedge or a contour shape to decrease the viscoelastic coating thickness and thereby minimize oscillations near the intersection of the viscoelastic coating with a rigid surface.

20. The method of claim 16, wherein the viscoelastic coating is combined with surface structures to enhance the stabilization of natural longitudinal and transverse wavelengths in the near-wall flow.

21. A viscoelastic coating designed according to the method of claim 16.

22. A viscoelastic coating designed according to the method of claim 17.

23. A viscoelastic coating designed according to the method of claim 18.

24. A viscoelastic coating designed according to the method of claim 19.

25. A viscoelastic coating designed according to the method of claim 20.

16. A method for reducing the friction drag of a body by providing a viscoelastic coating, said method comprising the following steps, performed in the order indicated:

a) determining characteristics of a turbulent boundary layer at a specified freestream velocity using boundary conditions for a rigid surface of an identical size and shape of the surface to be coated, said characteristics including viscous sublayer thickness, boundary layer thickness, phase speed and frequency corresponding to maximum-energy-carrying disturbances, mean velocity profiles, Reynolds stress distributions, wall shear stress distribution and friction drag;

b) selecting material properties of the viscoelastic coating including a density, a complex shear modulus and a coating thickness that, when excited by a forcing function substantially identical to that produced by the turbulent boundary layer as determined in step a), will provide the maximum flux of energy into the viscoelastic coating without producing surface excitation amplitudes that exceed the viscous sublayer thickness;

c) determining characteristics of the turbulent boundary layer over the viscoelastic coating and at the specified freestream velocity using oscillation amplitudes and an energy flux corresponding to the material properties of the viscoelastic coating selected in step b), including mean velocity profiles, Reynolds stress distributions, wall shear stress distribution, and friction drag,

d) determining a percentage reduction in friction drag as a ratio of the value of the friction drag without the viscoelastic coating as determined in step a) minus the friction drag with the viscoelastic coating as determined in step c) divided by the friction drag without the viscoelastic coating as determined in step a), and thereby quantifying the performance potential of the coating design; and

e) composing the viscoelastic coating from a material or a combination of materials as selected in all of the steps a) through d).