J Grandcoin

A Boukamel
email: boukamel@ehtp.ac.ma

S Lejeunes
email: lejeunes@lma.cnrs-mrs.fr

A micro-mechanically based continuum damage model for fatigue life prediction of filled rubbers

Keywords: Fatigue, Damage, Rubber, Finite Element

. A Finite Element implementation is proposed and some numerical examples are provided.

Introduction

Elastomeric materials are widely used by many industries: automotive, aeronautic... In many applications, elastomeric parts are closely linked to security and require specific qualification tests which impact product costs and safety of end-users. Therefore a perfect knowledge of the mechanical behavior and robust constitutive models are needed. Modern elastomeric components are composed of more and more complex rubber systems including different fillers and rubber bases (natural or synthetic or both). In many applications, these materials are subjected to cyclic loadings and are designed to damp energy. For these materials, it is mainly observed a damage initiation at a micro level by filler/matrix decohesion or by cavitation growth inside matrix. These phenomena give rise to a propagation at a macro level by the help of coalescence of voids or by crack propagation leading at the end to the failure of part. Identifying the physical phenomenon that takes place during high cycle fatigue test is a very hard task because multi-physics and multi-scale effects are occurring. For instance, temperature, crystallization, stress triaxiality, stress ratio, hydrostatic pressure, chemical aging are playing important roles on fatigue.

In the literature, one can distinguish at least three approach to model the damage due to cyclic loadings for rubber materials:

• The first one relates to prediction of crack growth. These approaches often assume that initial defects are pre-existing in the material but no information on defects localization or size is required. Starting from experimental curves that plot the value of a predicting quantity upon the log of number of cycles (so called Whöler curves or Haigh's diagrams) many predictor or criteria have been proposed. These criteria can be based on a mechanical quantity that is not directly related to local mechanisms caused by damage: maximum principal strain (e.g. [START_REF] Cadwell | Dynamic fatigue life of rubber[END_REF]), strain energy density (e.g. [START_REF] Beatty | Fatigue of rubber[END_REF]), strain based phenomenological model (e.g. [START_REF] Robisson | Comportement d'un élastomère sbr chargé de silice et d'un pu cellulaire: Prévision de la durée de vie en fatigue[END_REF]), stress based with Dang Van sphere (e.g. [START_REF] Brunac | On the heuristic extension of haigh's diagram for the fatigue of elastomers to arbitrary loadings[END_REF]), dissipated energy density (e.g. [START_REF] Lacroix | A local criterion for fatigue crack initiation on chloroprene rubber: approach in dissipation[END_REF]; [START_REF] Poisson | Biaxial fatigue behavior of a polychloroprene rubber[END_REF]). Other approaches try to embed more information on real local-mechanisms or idealized ones. [START_REF] Saintier | Crack initiation and propagation under multiaxial fatigue in a natural rubber[END_REF] have proposed a multi-axial criteria based on principal stresses with a reinforcing term linked to crystallization. In [START_REF] Mars | Cracking energy density as a predictor of fatigue life under multiaxial conditions[END_REF] a new criteria has been introduced, an incremental strain energy, which is defined on an unknown material plane that represents the cracking plane, is computed (see also [START_REF] Zine | Rubber fatigue life under multiaxial loading: Numerical and experimental investigations[END_REF]). [START_REF] Verron | A multiaxial criterion for crack nucleation in rubber[END_REF] have proposed a promising approach using a criteria based on the Eshelby tensor (configurational mechanics). This criteria shows very interesting results in multi-axial conditions (see also [START_REF] Verron | Definition of a new predictor for multiaxial fatigue crack nucleation in rubber[END_REF]; Andriyana et al. (2002)). Interesting comparisons of criteria can be found in [START_REF] Poisson | Biaxial fatigue behavior of a polychloroprene rubber[END_REF] for a polychloroprene rubber in multi-axial conditions and in [START_REF] Previati | Crack propagations in pneumatic tires: Continuum mechanics and fracture mechanics approaches[END_REF] for an application on truck tires. A more complete review of criteria can be found in [START_REF] Mars | Multiaxial fatigue of rubber[END_REF]; [START_REF] Mars | A literature survey on fatigue analysis approaches for rubber[END_REF].

• The second one, consists on using fracture mechanics to predict crack initiation and propagation.

One can find, for instance, power laws that describe the evolution of the crack length upon cycles as a power function of the energy release rate, see [START_REF] Lake | The mechanical fatigue limit for rubber[END_REF] for experimental results.

Fracture mechanics based models require a knowledge of initial defects which is very restrictive for rubber because these defects are very difficult to identify.

• The last approach is based on continuum damage mechanics (CDM) and consists in the introduction of, at least, one internal variable in the material behavior that represents the normalized quantity of damage (cracks, voids, etc) in the material. This approach has been firstly used for elastomeric materials to described short time history effect like the Mullins effect (see [START_REF] Chagnon | On the relevance of continuum damage mechanics as applied to the mullins effect in elastomers[END_REF]).

For high cycle fatigue, the ideas of [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF] can be applied to rubber. A cumulative rule for damage is assumed (e.g. linear) and the damage evolution is integrated cycle by cycle. Therefore, the damage evolution can be directly related to the number of cycles. This approach has been used in [START_REF] Wang | A damage model for the fatigue life of elastomeric materials[END_REF]; [START_REF] Ayoub | A continuum damage model for the high-cycle fatigue life prediction of styrene-butadiene rubber under multiaxial loading[END_REF][START_REF] Ayoub | Fatigue life prediction of rubber-like materials under multiaxial loading using a continuum damage mechanics approach: Effects of two-blocks loading and r ratio[END_REF]. To the authors knowledge only hyperelasticity has been considered in previous works, however the mechanical dissipation plays a fundamental role on the fatigue life in rubber materials.

In the present work a silicone rubber filled with silica is considered. This material is used in aeronautic structures to damp a mechanical energy due to vibrations. In industrial applications, rubber is mainly loaded in simple shear and the prediction of a stiffness loss is critical for the safety of the structure. It is proposed to use the CDM approach to predict this stiffness loss and the fatigue life. Starting from micromechanical observations and from an original experimental campaign of fatigue tests on simple shear and tension specimens, we propose a two scale model: the first scale aims at representing the voids or cracks initiation near the agglomerate of fillers and the second scale represents the voids or cracks growth in the rubber matrix. The main originality of this paper is that the proposed driving force for damage is defined from a returnable energy rate which is integrated over a cycle. This energy represents the difference between the given energy and the dissipated one for a stabilized cycle. The effect of stress multi-axiality or the effect of the mean stress are embedded in this energetic damage rule.

The constitutive behavior developed in this paper is based on the statistical framework of [START_REF] Martinez | Statistical approach for a hyper-visco-plastic model for filled rubber: Experimental characterization and numerical modeling[END_REF]. In this framework, the Cauchy stress is viewed as the sum of a matrix part (hyperelastic) and a continuous contribution due to agglomerates (visco-hyperelastic). A statistical variable and a probability function are introduced to represent a continuous number of relaxation mechanisms. The statistical variable is related to the activation energy of each dissipation mechanisms (related to each population of agglomerates). The main advantage of this approach lies in its ability to cover a large range of strain rates and amplitudes with a small number of material parameters. The connection between the damage model and the constitutive model is straight-forward, the micro scale of damage is taken into account in the probability function as damage in agglomerates tends to decrease the relaxation mechanisms. The macro scale is classically related to the cracks grow in the matrix.

A classical variational formalism is adopted and a finite element implementation is derived. The mechanical flow rules are integrated following the proposal of [START_REF] Lejeunes | Finite element implementation of nearly-incompressible rheological models based on multiplicative decomposition[END_REF]. Damage is only post-treated by computing the returnable energy rate over one stabilized cycle.

The paper is organized as follows. In a first part, the micro-mechanical observations and the main results of the experimental characterization are presented and discussed. In a second part, a new statistical constitutive model is introduced to represent the mechanical behavior. The CDM approach is then detailed in a third section. The finite elements implementation is briefly exposed and some numerical results are shown in the last section. To identify micro-mechanisms that take place during cyclic fatigue test, we made some SEM observations at the "Centre des Materiaux" of the "Ecoles des Mines de Paris". Figure 1(a) shows a large agglomerate of (a) matrix/filler decohesion at 100% of uniaxial-tension: the tension direction is in the same direction as the one of loading (horizontal direction for this image) (b) matrix/filler decohesion at 0% of uniaxial-tension: cavitation at the poles of the agglomerate Figure 2: SEM observations of silica agglomerates in a silicone rubber after fatigue: In-Situ tension test. Due to the metalic agent to deposit on the surface of the specimen prior to In-Situ testing, the cracks observed can not be directly related to the rubber matrix.

Experimental results

micro-mechanical observations

silica inside the matrix for a virgin sample, at this scale (20µm) the matrix is dense, no voids are observed. The figure 1(b) shows a smaller agglomerate for a sample which has been cut inside a uni-axial tension specimen. This specimen has been previously submitted to 100000 cycles of fatigue (frequency is 1Hz and the strain amplitude is 150% ± 50%). A decohesion of the filler-matrix interface can be observed. On figure 1(c) and 1(d), one can observed the initiation of a micro-crack (5 µm length on figure 1(c)) which is normal to the direction of the fatigue loading. This micro-crack contains a filler at its center and it seem quite natural to postulate that the filler can be viewed as a defect that originated the crack, however we have no proof that the crack was not previously here. The SEM also provide the availability to do In-situ observations (while applying a mechanical loading). Figures 2 show a result of such an experiment. The micro sample, which has been cut in a previously fatigued specimen, is submitted to an uni-axial tension test (the loading direction is as the one of fatigue). A filler/matrix decohesion is seen, in the same direction as the one of loading. When the load is released a void remains at the poles of the agglomerate (even after a long time to eliminate viscous effects).

From these observations we can conclude that filler/matrix interaction plays a fundamental role on the fatigue behavior at the micro-scale for this material. Decohesion seems to be the principal micro-mechanism that can lead to the formation of micro-cracks. For uni-axial loading and for a non-cristallizable matrix, these micro-cracks propagate in the matrix in a direction which is normal to the load direction without bifurcations. Furthermore, this decohesion phenomena leads to a volume variation because voids remain in the unloading state. However, transposing these results at the macro-scale is not obvious. Therefore we prefer to neglect it in the following as we believe that at the macro-scale damage due to fatigue is affecting more the deviatoric part than the spherical one. To the authors opinion the investigation of volume variation is a very hard topic as concurrent phenomena are taking place in rubber that affect the volume variation: cavitation, crystallization, network reorganization or stress softening, thermal dilatation, etc.

fatigue characterization

The aim of the experimental campaign is to characterize the mechanical behavior after a given number of cycles. Usually, damage is characterized through a loss of rigidity. However, rubber materials exhibit other softening effects like the Mullins effect or thermal softening induced by self-heating for high dissipative materials.

In fact, the increase of temperature during the fatigue test implies a softening of the material. If the acquisitions are done during fatigue tests, the thermal softening can hide the damage effect which also leads to a softening. To clearly identify the evolution of damage with a remove of the thermal and the Mullins softening, we have defined an original experimental protocol. The characterization of the behavior at a given number of fatigue cycles N requires the following steps:

1. 20 cycles at a strain amplitude 10% higher than the maximum strain amplitude during fatigue to remove Mullins effect, 2. N fatigue cycles with a regulated temperature (only the ambient temperature inside the climatic chamber is controlled), 3. two rest days (the sample is stress free) to remove thermal and long time viscous effects, 4. mechanical characterization with: (i) a quasi-static test, (ii) triangular cyclic tests at various strain rates and various strain amplitudes and (iii) relaxation tests.

All these tests have been performed on a hydraulic machine which can be equipped with an adiabatic chamber in which internal temperature can be regulated. For shear specimens, displacement is controlled and load is recorded. For tension specimen depending on the load case, the control can be done on displacement or on load for the unloading phase. Therefore it can be avoided buckling due to relaxation. The signal is sinusoidal. In the following, for each fatigue test it is provided the frequency of the signal, the maximum strain amplitude, and the mean strain if it is not zero. Tension specimens are dumbbell shaped of type H2 with a square rectangular cross-section of 5x2 mm 2 . For shear tests, it is used specific double shear specimens. These specimen are composed of two blocks of rubber with a central rigid aluminum layer and two external rigid aluminum layers. The rubber blocks are 96 mm long, 5 mm high and 20 mm large. Furthermore the end sections have a meniscus shape to ensure that cracks do not initiate at the interface aluminum/rubber. Each specimen is submitted further time to the steps 1 to 4 with increasing value of N (typically N=1000, 10000, 100000,...)

The previous procedure is time consuming and only a repeatability of two was achieved for simple shear specimens and three for tension specimens. Figures 3, show an example of experimental results obtained from the step 4 (characterization tests). In figure 3(a) and 3(d) we clearly see a loss of stiffness and a decrease of the dissipated energy when N (the number of fatigue cycles) increases. More surprisingly, we often observe during the first cycles of fatigue in simple shear tests an increase of the stiffness (see figure 3(b) and 3(c)) previously to damage phenomena. This observation will be discussed below.

To proceed further, it is defined global characteristics (stiffness and dissipation) from the figure 4. Figures 5 and6 show a comparison of different fatigue configurations with tension and simple shear specimens. Each point is obtained by following all the steps previously described. For each quasi-static characterization tests after fatigue, we determine a global stiffness and a global dissipation on a stabilized cycle. The stiffness and the dissipation values are normalized and it is observed the following behavior: in a first step the global stiffness and the dissipation are increasing followed by a second step in which both stiffness and dissipation are rapidly decreasing. Furthermore this effect is more pronounced for shear specimens. We propose the following scenario to explain these observations:

Step 1 During the first cycles, fillers-matrix bonds begin to break and the filler-matrix decohesion is starting. In the meantime, new chemical bonds are formed in the matrix due to the thermal aging that results from the self-heating of the material. These effects are in competition but thermal aging takes the lead and a global stiffening is observed.

Step 2 For a large number of cycles, damage is growing and propagates to the matrix. Thermal aging is less sensitive due to the stabilization of the self-heating phenomena. A global decrease of the stiffness is observed.

The previous scenario also furnishes an explanation to the fact that shear specimens are more sensitive to thermal aging: in simple shear damage is more localized and less homogeneous than for tension and self-heating phenomena is stronger. Due to the competition between damage and thermal aging it is a very hard task to determine when initiation of damage in the matrix takes place. To proceed further we propose to consider the intersection of the two curves as an important characteristic number of cycles for which it can be clearly said that the mechanical behavior is affected by damage. Different fatigue configurations are summed up in the table 1. The loss of stiffness (and dissipation) is strongly dependent on the strain amplitude and it seems also sensitive to the mean strain: increasing R-ratio leads to higher number of cycles. This effect may be related to previous works: e.g. [START_REF] Legorju-Jago | Fatigue initiation and propagation in natural and synthtic rubbers[END_REF]; [START_REF] Saintier | Cyclic loadings and crystallization of natural rubber: An explanation of fatigue crack propagation reinforcement under a positive loading ratio[END_REF] show that life-time increases for positive R-ratio on diabolo specimens. They show that cristallization plays an important role on it but

Constitutive model

Mechanical behavior

In this paper, we extend a previously developed approach for the modeling of the dynamical behavior of filled rubbers (see [START_REF] Martinez | Statistical approach for a hyper-visco-plastic model for filled rubber: Experimental characterization and numerical modeling[END_REF]). This approach is micro-mechanically motivated and is based on a statistical representation of agglomerates and aggregates of fillers. Each population of aggregate or agglomerate is associated to a couple of thermodynamical potentials: a specific free energy and a pseudopotential of dissipation. A random variable ω, is introduced and relates to the energy of activation of a micro-mechanism of relaxation. It is also introduced P (ω), a probability density function that describes the range of energies of activation that is covered by the material. Therefore it is possible to take into account of a large number of relaxation mechanisms without increasing the number of material parameters.

It is assumed that the global deformation of the matrix and of a representative volume element of filler aggregates1 are the same (affine model). The global framework can be summarized by the figure 7. In this work, the matrix is described by a visco-elastic model of Poynting-Thomson. This viscosity is introduced to represent long time effects. Each population of agglomerate is described by an elasto-visco-plastic model. Following Flory [START_REF] Flory | Thermodynamic relations for highly elastic materials[END_REF]), the deformation gradient F is split into an isochoric part: F and a volumetric part: J 1/3 1, where J = det F and 1 stands for the identity tensor (Cartesian metric is used). As no viscous effects are generally observed for volumetric deformation in rubber, the isochoric deformation gradient is split as follows:

F = F e Fi = fe (ω) fi (ω) (1)
where F e and Fi are elastic and inelastic isochoric deformation gradients in the matrix and fe (ω) and fi (ω) are elastic and inelastic isochoric deformation gradient in the agglomerates. The themodynamical potentials

(free energy and pseudo potential of dissipation) are therefore functions of elastic and inelastic gradients but also depend on ω as follows:

     ψ = ψ matrix ( Bi , Be ) + ∞ 0
ψfillers ω, be (ω) P (ω)dω + ψ vol (J)

ϕ = ϕ matrix ( Do i ) + ∞ 0 φfillers ω, do i (ω) P (ω)dω (2)
where Be = F e F e T , be (ω) = fe (ω) fe (ω) T , Bi = Fi FT i are elastic and inelastic left Cauchy-Green tensors, Do i and do i (ω) are objective rates of the inelastic deformations defined from:

Do i = R e Di R e T = 1 2 R e ( Ḟi F-1 i + F-T i Ḟi T )R e T ( 3
)
do i (ω) = r e (ω) di (ω)r e T (ω) = 1 2 r e (ω)( ḟi (ω) f -1 i (ω) + f -T i (ω) ḟi T (ω))r e T (ω) (4)
where R e and r e (ω) are coming from the polar decomposition of F e and fe (ω).

Neglecting thermal effect, the dissipation inequality (or Clausius Duhem inequality) in eulerian configuration can be written as:

φ = σ : D -ρ ψ ≥ 0 (5)
The variation of the free energy, eq. ( 2), is given by:

ψ = ∂ψ matrix ∂ Be : Ḃe + ∂ψ matrix ∂ Bi : Ḃi + ∞ 0 ∂ ψfillers ω, be (ω) ∂ be (ω) : ḃe (ω)P (ω)dω + ∂ψ vol ∂J J (6)
The time variation of elastic and inelastic Cauchy Green tensors and of volume variation are given by:

J = J (1 : D) (7)
Ḃi = Li Bi + Bi LT i (8) Ḃe = L Be + Be L T -2 Ve Do i Ve - 2 3 (1 : D) Be (9) ḃe (ω) = L be (ω) + be (ω)L T -2v e (ω) do i (ω)v e (ω) - 2 3 (1 : D) be (ω) (10)
where Ve and ve (ω) are the pure deformations coming from the polar decomposition of F e and fe (ω). Using equations ( 10), ( 9), ( 8), ( 7) in ( 6) and reporting the result in (5), it is obtained the following expression of the dissipation:

φ =   σ -ρ 0 J -1 2 Be ∂ψ matrix ∂ Be D -ρ 0 J -1 ∞ 0 2 be (ω) ∂ ψfillers ∂ be (ω) D P (ω)dω -ρ 0 ∂ψ vol ∂J 1   : D + 2ρ 0 J -1 Ve ∂ψ matrix ∂ Be Ve -2ρ 0 J -1 V-1 e F e Bi ∂ψ matrix ∂ Bi F e T V-1 e : Do i + 2ρ 0 J -1 ∞ 0 ve (ω) ∂ψ f illers ∂ be ve (ω) : do i (ω)P (ω)dω ≥ 0 (11)
From equation ( 11), it is deduced the following constitutive equation for stress:

σ = σmatrix 2ρ 0 J -1 Be ∂ψ matrix ∂ Be D + σ f illers ∞ 0 2ρ 0 J -1 be (ω) ∂ ψfillers ∂ be (ω) D σfiller (ω) P (ω)dω + σ vol ρ 0 ∂ψ vol ∂J 1 (12)
It remains two part in the dissipation which are assumed to be independently positives. A quadratic pseudopotential is chosen for the matrix:

ϕ matrix ( Do i ) = η 2 ( Do i : Do i ) (13)
Where η is a viscosity parameter of the matrix. Therefore, using the concept of normality, the thermodynamical force that drives inelastic effects in the matrix is defined from the variation of the pseudo-potential.

Assuming the isotropy of the free energy function, the first complementary equation is obtained:

Do i = 2ρ 0 Jη Be ∂ψ matrix ∂ Be -2 V-1 e ψ matrix,1 B -ψ matrix,2 Be B-1 Be V-1 e D (14)
In the previous equation, ψ matrix,1 and ψ matrix,2 stands for the variations of the isotropic free energy upon the inelastic invariants:

ψ matrix,1 = ∂ψ matrix ∂I 1 ( Bi ) ψ matrix,2 = ∂ψ matrix ∂I 2 ( Bi ) ( 15
)
Concerning the second part in the residue of dissipation, as P (ω) is a positive function the following condition is sufficient to fulfill the Clausius-Duhem inequality:

ve (ω) ∂ψ f illers ∂ be ve (ω) : do i (ω) ≥ 0 ∀ω ∈ [0, ∞[ (16)
We define an elastic domain:

E(ω) = { σfiller |f ( σfiller , ω) ≤ 0} ( 17
)
where f is a yield stress function for which the yield depends on ω. We assume the following form for the dual2 pseudo-potential of dissipation:

ϕ * = f ( σfiller , ω) 2 2β(ω) with f ( σfiller , ω) = σfiller -χ(ω) (18)
where β(ω) is a viscosity function, < . > are the Mac-Cauley brackets3 , χ(ω) is the yield parameter. Using the normality principle, we obtain the following expression:

do i (ω) = ∂ϕ * ∂ σfiller = f ( σfiller ) β(ω) σfiller σfiller (19)
To proceed further, we make the following choices for the free energies:

   ρ 0 ψ matrix = C 1 I 1 ( Be ) -3 + C 2 Ln I 2 ( Be ) 3 + A I 1 ( Bi ) -3 ρ 0 ψfillers ω, be (ω) = a(ω) I 1 ( be (ω)) -3 (20)
where C 1 , C 2 , A are material parameters. Concerning the statistical functions, we choose a Gaussian form for the density of probability, a linear form for the yield parameter and the viscosity of the fillers, the shear modulus of the fillers follows an exponential decay: To resume, the constitutive model is defined from the following set of constitutive equations and flow rules:

             P (ω) = 1 P 0 exp(-( ω Ω ) 2 ) with P 0 = ∞ 0 exp - ω Ω 2 dω χ(ω) = χ ω ω0 , a(ω) = a 0 exp -ω ω0 , β(ω) = β ω ω0 (
σ = σ matrix + ∞ 0 σfiller P (ω)dω + p1 σ matrix = 2J -1 C 1 Be + C 2 I 2 ( Be ) (I 1 ( Be )1 -Be ) D σfiller = 2J -1 a(ω) be (ω) D p = ρ 0 ∂ψ vol ∂J Ḃe = L Be + Be L T - 2 3 (1 : L) Be - 2 η σ matrix Be + 4A Jη B - 1 3 ( B : B-1 e ) Be ḃe (ω) = L be (ω) + be (ω)L T - 2 3 (1 : L) be (ω) -2 f ( σfiller ) β(ω) σfiller σfiller be (ω) (22)
There are 9 material parameters to identify, 7 are deterministic (C 1 , C 2 , A, η, a 0 , β and χ) and 2 are statistical (Ω and ω 0 ).

Continuum fatigue damage model

We start from the following hypothesis: (1) the damage is isotropic, (2) the damage is accumulated in a linear manner and (3) the damage is composed by two contributions: (i) a microscopic one, that corresponds to the fillers/matrix de-cohesion observed with the SEM and (ii) a macroscopic one which is activated upon reaching a yield damage that corresponds to cracks propagation in the matrix leading to a stiffness loss.

In the constitutive model the statistical parameter Ω controls the probability to activate a micromechanism of dissipation in the aggregates. A larger value of Ω leads to a larger standard deviation. Therefore when matrix/fillers decohesion occurs the stress contribution of agglomerates decreases which can be represented by a smaller standard deviation. The micro damage d is then defined from the following relation:

Ω = (1 -d) Ω ( 23
)
where Ω is the parameter of a virgin material. The macro damage D which concern only the matrix is classically defined from:

σ matrix = (1 -D) σmatrix ( 24
)
where σmatrix is the effective stress (which correspond to a virgin material). The figure 8 shows a schematic representation of the supposed damage scenario. The internal variable rates ḋ and Ḋ can be derived from a pseudo-potential of dissipation defined from the thermodynamical forces associated to damage, however in the following we consider the damage increment over one cycle such as: In this approach we postulate that the accumulation of damage is related to the difference between the power of internal forces and the dissipated energy integrated over a cycle. This energetic quantity is referred as a returnable cyclic energy. We define R matrix and R f illers , respectively the returnable cyclic energy of matrix and agglomerates such as:

       R matrix = cycle ρ ψmatrix dt = cycle σ matrix : D -φ matrix dt R f illers = ∞ 0 cycle ρ ˙ ψ f illers (ω) dt P(ω)dω = cycle σ f illers : D -φ f illers dt (26)
where < x > denotes the positive part of x. Starting from these energies, we state a power law evolution for the micro damage over a cycle:

δd δN = d α 1 -α R f illers R 0 f illers β with α < 1 (27)
where α, β and R 0 f illers are material parameters. For the macro damage over a cycle we also assume a power law evolution but with a yield function depending on the level of micro-damage. Therefore, we assume that a macro damage could not initiate if a micro-damage do not have reached the yield limit d 0 :

δD δN = H(d -d 0 ) D ζ 1 -ζ R matrix R 0 matrix ν with ζ < 1 (28)
with H(x) a Heaviside function and ζ, ν, d 0 , R 0 matrix material parameters. The integration of the equation ( 27) and (28) over cycles leads to the following:

d =   N R f illers R 0 f illers β   1 1-α (29) D = H(d -d 0 ) N R matrix R 0 matrix ν 1 1-ζ (30)
From the previous equations, it can be computed two maximum number of cycles N micro limit and N macro limit , which are defined from d = 1 and D = 1:

                   N micro limit = R f illers R 0 f illers -β N macro limit = R matrix R 0 matrix -ν if N macro limit > d 1-α 0 R f illers R 0 f illers -β N macro limit = ∞ otherwise (31)
The end of life is reached when D = 1 which corresponds to the occurence of a macro crack in the matrix. The propagation of these macro cracks must be described by fracture mechanics and are not considered in this paper.

Parameters identification

Undamaged behavior

We use the same strategy as proposed in [START_REF] Martinez | Statistical approach for a hyper-visco-plastic model for filled rubber: Experimental characterization and numerical modeling[END_REF]. This strategy consists in three steps: (1) Quasi-static experimental tests (tension and simple shear) allow to determine C 1 , C 2 and to estimate A, a 0 , χ, Ω and ω 0 ; (2) Relaxation tests (tension and simple shear) allow to identify A, H and ω 0 , and to correct a 0 , χ, η, Ω and (3) triangular cyclic tests, at different rates (simple shear) allow to correct the prediction of a 0 , χ, η and Ω.

Each step of identification is realized with an heuristic algorithm of minimization of the sum of the squared distances between the experimental data and the analytical or semi-analytical responses. The material parameters obtained are presented in table 2.

Figures 9 show the responses of the identified model compared to experimental tests. Some tests have been used in the identification strategy, this is the case for figures 9(a), 9(b) and 9(c). Others are used for validation purpose, this is the case of figure 9(d). Comparatively to the previous version of the statistical approach that has been presented in [START_REF] Martinez | Statistical approach for a hyper-visco-plastic model for filled rubber: Experimental characterization and numerical modeling[END_REF], the present model exhibits a better agreement to the experimental tests.

Damage parameters

To identify the damage parameters, only uni-axial tension tests are used because we assume that the central zone of the specimen sees a uniform damage which is clearly not the case of the simple shear specimen.

13

C 1 C 2 A η a 0 β χ Ω ω 0 1.
28 Mpa -2.67 Mpa 2.95 Mpa 50.13 Mpa.s 5.12 Mpa 0.5 Mpa.s 1.69 Mpa 3.11e -2 0.40 Therefore, for each fatigue test in tension it is identified a set of material parameters as explained in the previous section. Using the definition of eqs. ( 23) and ( 24), it is computed the value of d, D, R matrix and R f illers for each configuration of fatigue. Taking logarithmic value of eqs. ( 27) and ( 28), one can separate the identification of α and A from the others parameters. Figures 10(a) and 11(a) show the evolution of the damage rate upon the current value of damage computed after each characterization test done after a fixed number of cycles (repeated for each configuration of fatigue). The parameters α and A can be identified from these curves. Figures 10(b) and 11(b) show the dependance of damage upon returnable cyclic energy. The parameters β, B, R 0 f illers and R 0 matrix can be computed from these curves. The last parameter, d 0 , is determined by computing a mean of d values for which D starts to increase. The parameters are given in table 3.

Finite elements simulations

Variational Formulation

To resolve the equilibrium equations in a quasi-static context with a nearly-incompressible constitutive behavior it has been chosen a perturbed Lagrangian weak form. This formulation is based on the introduction of a Lagrange multiplier, p which stands for the hydrostatic pressure. Therefore, the equilibrium solution is defined by the couple (u, p) which have to cancel the following integral form for all the test functions (δu, δp):

       Ω0 (Π matrix + Π f illers + pJF -T ) : ▽δudΩ - Ω0 δu • f vol dΩ - δΩ0 f δu • f surf dS = 0 Ω0 g(J) - p K δpdΩ = 0 (32)
where Ω 0 is the initial domain (in the reference configuration), δΩ 0f is the initial boundary where current surface forces f surf are applied, f vol are the current volumetric forces. Tensors Π matrix and Π f illers stand for the first Piola-Kirchof stress respectively in the matrix and in the fillers, they are defined from the following relations:

Π matrix = Jσ matrix F -T and Π f illers = Jσ f illers F -T (33)
The parameter K is the compressibility modulus and its inverse plays a role of a perturbation parameter. The compressibility function g(J), is defined from the following relation: 34) in this paper, the classical linear compressibility law is used leading to the following free energy potential:

g(J) = ρ 0 1 K ∂ψ vol ∂J (
ρ 0 ψ vol = K 2 (J -1) 2 (35)

Numerical Implementation

The integration of the constitutive model defined by eqs ( 22), requires to compute the statistical integral form and the integration of the local flow rules. The statistical integration is realized by a classical trapezoidal scheme such as:

σ f illers = n-1 k=1 1 2 ( σfillers (ω k )P (ω k ) + σfillers (ω k+1 )P (ω k+1 )) (ω k+1 -ω k ) (36)
where ω k is defined from the number of integration points n such as: ω k = k(ω e )/n, ω e is an activation energy limit associated with a probability of 0.01P 0 :

ω e = -ln(0.01)Ω 2 (37)
Therefore, the computation of σ f illers requires the resolution of n flow rules. Each flow rule, is integrated with a specific backward Euler scheme which is detailed in [START_REF] Lejeunes | Finite element implementation of nearly-incompressible rheological models based on multiplicative decomposition[END_REF]. In practice, these computations are easy to parallelize as the flow rules are not dependent from each other. Therefore, it has been chosen the following strategy: at each gauss point the sum defined by eq. ( 36) is parallelized using an OpenMP framwork. Each cpu core solves a flow rule with its own copy of b e (ω). The computing time is directly related to the stiffness of the differential system of equations obtained for each flow rule. Therefore this strategy is not optimal but it is very easy to implement. The computation of a local state of damage requires the evaluation of the returnable cyclic energy R matrix and R f illers , as for the statistical integration, a trapezoidal integration rule is used to compute the integrals of equation ( 26).

The space discretization is realized with classical stable mixed elements with a linear interpolation for the pressure and a quadratic one for the kinematic field.

Simulation of the simple shear test

As mentioned previously, the simple shear test leads to an inhomogeneous state of damage, therefore it is proposed to use the results of this test as a validation purpose. As mentioned previously, the shear specimens used in this work, exhibit a specific shape with a meniscus on each lateral faces to enforce an initiation of damage in rubber and not in the interface rubber/aluminum (see figure 12(a)). Experimentally it was observed a macroscopic crack initiation near the center of the meniscus which propagates with an approximate angle of 45 o . In some pathological tests the crack rapidly attains the aluminum layer and a decohesion of the aluminum/rubber interface is observed (see figure 12(b)).

For the finite elements simulation we assumed plane strain conditions, the aluminum is supposed to be perfectly rigid and not modelized. The mesh comprises 1183 nodes and 265 quadratic elements. To ensure that the numerical response corresponds to a stabilized response it is done five cycles of shear loading (transversal displacements with no vertical displacements are imposed on the top edge). The returnable cyclic energies are computed on the last cycle and the maximum number of cycles (lifetime) is deduced.

The figure 13 shows the evolution of the two damage variables, d and D in the shear specimen for a fatigue configuration of ±50% of shear strain at 6Hz. It can be observed first that the maximum of macro damage is located near the center of the meniscus which is in good accordance to experimental observations. Furthermore, it can be also seen that the macro damage seems more inhomogeneous than the micro damage. Figure 14 shows a comparison of the global response of the numerical shear specimen and the experimental data. The stiffness loss seems well predicted by the model. It can also be computed a predicted life time (N macro limit ) which is given by the minimum number of cycles for which D = 1 in a sufficiently large zone of the specimen (here it has been considered that 90% of the specimen area was a sufficiently large zone.). The table 4, presents some results of predicted life time for different fatigue configurations. As already seen experimentally it is observed a non linear dependence upon the solicitation amplitude and a weak influence of the frequency, at least in the experimental range. Unfortunately, these results cannot be confronted to experimental one as we do not have made enough tests in simple shear to be representative (the aluminum/rubber interface has broken suddenly in many tests).

Concluding remarks

The use of continuum damage mechanics for rubber seems promising for at least two reasons. First, the experimental campaign for such an approach is less time consuming than approaches that are based on Whöler curves or Haigh's diagram. Furthermore, interlaying characterization tests in a fatigue test allows to separate thermal softening from damage softening. Second, the model is not restricted to only one damage mechanism or one physical criteria to predict the fatigue behavior. In the present paper, it has been considered two damage mechanisms which are micro-mechanically motivated: filler agglomerates decohesion and matrix voids or cracks growth. Therefore, these models can propose a link between micro-mechanical observations and fatigue behavior.

The silicon rubber that is considered in this paper exhibits a complex behavior. We have observed a thermal aging previously to damage softening. We assume that this effect is due to a chemical activity (cross-liking) during the fatigue test. However, a systematic chemo-physical characterization is needed to study this effect.

The proposed model is based on a micro-motivated approach with a statistical representation of agglomerates populations. It allows a continuous representation of the characteristic times of viscosity. Therefore, the introduction of a micro damage variable that affects fillers/matrix interface behavior is straight forward: only one statistical variable has to be considered. The main originality of this model resides in the introduction of the returnable cyclic energy to build damage evolution laws. This energy takes into account of dissipation and free energy. In this paper, it is only considered mechanical dissipation as thermal effects has been neglected. However, the proposed criteria can be extended to a thermo-mechanical case.

The finite element implementation of the proposed model has allowed to confront some numerical results in simple shear to experimental ones for a model which has been identified in tension. First results are encouraging but more experimental tests have to be done. In a future work, we will looking at diabolo specimens to do tension/torsion or compression/torsion tests.

  (a) Virgin material (large agglomerate of silica) (b) After fatigue: void near a filler (probably due to matrix/filler decohesion) (c) After fatigue: micro-crack normal to the load direction (d) After fatigue: zoom on micro-crack with ligaments between filler and matrix.

Figure 1 :

 1 Figure 1: SEM observations of Silica agglomerates in a silicone rubber.

  Triangular cyclic tension tests after N cycles of fatigue (the fatigue configuration is 0.5Hz and 150% of mean strain and ±150% of dynamic strain). Relaxation tests in tension after N cycles of fatigue (the fatigue configuration is 0.5Hz and 25% of mean strain and ±25% of dynamic strain). Triangular cyclic simple shear tests with progressive amplitude, after N cycles of fatigue (the fatigue configuration is 15Hz and 12.5% of mean strain and ±25% of dynamic strain). Triangular cyclic simple shear tests with progressive amplitude, after N cycles of fatigue (the fatigue configuration is 15Hz and 12.5% of mean strain and ±25% of dynamic strain).

Figure 3 :

 3 Figure 3: Typical results obtained from the characterization tests after fatigue (step 4 of the experimental campaign)

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Global characteristics: global stiffness (slope of the secant of the hysteresis) and global dissipation (area of the hysteresis)

Figure 8 :

 8 Figure 8: Supposed damage senario

Figure 9 :

 9 Figure 9: Comparison of the identified model and some experiments in shear and tension.

  Damage rate evolution for different fatigue tests at 0.5Hz ( : 25±25% of strain, + : 50±50% of strain, × : 100±100% of strain, : after 150 ± 150% of strain).

  Figure 10: Micro damage identification

  Figure 11: Macro damage identification

  Figure 12: Shear specimen

Figure 14 :

 14 Figure 14: Evolution of the hysteresis upon the number of fatigue cycles for the shear specimen. Sinusoidal shear fatigue test at γ = 0.5 and f = 6Hz. Frequency (Hz) Shear amplitude (%) N macro limit (number of cycles) 6 12.5 2044793 6 25 931663 6 50 400564 15 25 931779 15 50 401318

  [START_REF] Poisson | Biaxial fatigue behavior of a polychloroprene rubber[END_REF] also observed it for non-cristallizable rubber. Frequency does not seems to play a major role but more tests are needed to study this effect.

		Maximum strain amplitude (with zero	Maximum strain amplitude (with 12.5% of
		mean strain)			mean strain)
		12, 5%	25%	50%	12, 5%	25%	50%
	6Hz	> 1000000	168000	30000	> 1000000	190000	[30000, 50000]
	15Hz > 1000000	148000	30000	> 1000000	191000	[30000, 50000]

Table 1: Number of cycles for which global stiffness and global dissipation start to decrease for the simple shear fatigue test global stiffness global dissipation

Table 2 :

 2 Material parameters

Table 3 :

 3 Damage parameters

Table 4 :

 4 Predicted life time (D = 1) for different fatigue configuration in simple shear

the RVE of aggregates is assumed to be composed of fillers and gum (occluded gum and bound rubber).

Dual pseudo-potential (which depend on force) can be obtained from a Legendre-Fenchel transformation of the primal pseudo-potential (which depend on flux)

Operator < . > is defined by < f >= (f + |f |)/2


Ashley P Willis
email: willis@ladhyx.polytechnique.fr

Yongyun Hwang
email: yongyun@ladhyx.polytechnique.fr

Carlo Cossu
email: carlo.cossu@imft.fr

Optimally amplified large-scale streaks and drag reduction in turbulent pipe flow

The optimal amplifications of small coherent perturbations within turbulent pipe flow are computed for Reynolds numbers up to one million. Three standard frameworks are considered: the optimal growth of an initial condition, the response to harmonic forcing and the Karhunen-Loève ͑proper orthogonal decomposition͒ analysis of the response to stochastic forcing. Similar to analyses of the turbulent plane channel flow and boundary layer, it is found that streaks elongated in the streamwise direction can be greatly amplified from quasistreamwise vortices, despite linear stability of the mean flow profile. The most responsive perturbations are streamwise uniform and, for sufficiently large Reynolds number, the most responsive azimuthal mode is of wave number m = 1. The response of this mode increases with the Reynolds number. A secondary peak, where m corresponds to azimuthal wavelengths + Ϸ 70-90 in wall units, also exists in the amplification of initial conditions and in premultiplied response curves for the forced problems. Direct numerical simulations at Re = 5300 confirm that the forcing of m = 1 , 2 and m = 4 optimal structures results in the large response of coherent large-scale streaks. For moderate amplitudes of the forcing, low-speed streaks become narrower and more energetic, whereas high-speed streaks become more spread. It is further shown that drag reduction can be achieved by forcing steady large-scale structures, as anticipated from earlier investigations. Here the energy balance is calculated. At Re= 5300 it is shown that, due to the small power required by the forcing of optimal structures, a net power saving of the order of 10% can be achieved following this approach, which could be relevant for practical applications.

I. INTRODUCTION

The nature of laminar-turbulent transition, and of the structure of turbulence in pipe flow, is a subject that has interested generations of fluid dynamicists since the pioneering work of Reynolds ͓1͔. Laminar pipe flow ͑the Hagen-Poiseuille flow solution͒ is linearly stable, but subcritical transition to turbulence is observed at Reynolds numbers as low as Ϸ2000. A key mechanism in the transition process is played by the lift-up effect ͓2,3͔ consisting of the strong amplification of the energy of streamwise streaks that evolve from streamwise vortices. These streaks are narrow elongated regions where the streamwise velocity is larger or smaller that the average at the same distance from the wall. The large amplification of the streaks by the lift-up mechanism is related to the strong non-normality of the linearized Navier-Stokes operator ͓4-6͔. Optimal perturbations, which maximize the energy amplifications, have been computed for most of the canonical laminar shear flows, where it is found that streamwise uniform structures are the most amplified. The maximum amplification of initial conditions, and responses to harmonic and stochastic forcing have been found to scale such as Re 2 , Re 4 , and Re 3 , respectively ͓6-9͔. For laminar pipe flow, the most amplified streamwise uniform modes are those of azimuthal wave number m = 1 ͓10-12͔, decreasing monotonically with increasing m, and the resulting streaks are significantly unstable when of sufficient amplitude ͓13͔.

Coherent streaky motions also exist in the turbulent pipe. The existence of streaks in the buffer layer with a mean spanwise spacing of Ϸ100 wall units has been well documented for the boundary layer ͓14-16͔ and for pipe flows ͓17,18͔. Coherent streaky structures exist also at much larger scales, of the order of several pipe radii, R. Recently it has been realized that these structures can be very long ͑at least Ϸ8R to 16R͒, and contribute to half of the total turbulent kinetic energy and Reynolds stress ͓19,20͔. Contrary to other canonical turbulent flows, for the case of the turbulent pipe relatively little is known on the spanwise ͑azimuthal͒ structure of these coherent large-scale structures.

While early computations of optimal growths and responses to stochastic forcing ͓21-23͔ have considered the turbulent mean flow as a base flow, the molecular viscosity was included in the equations for the coherent perturbations. An externally imposed eddy turnover time had to be included in the optimization process. More recent analyses of the turbulent plane channel flow ͓24-26͔, boundary layer ͓27͔ and Couette flow ͓28͔ have included the effect of the eddy viscosity in the linearized equations, in the spirit of Ref. ͓29͔. In this approach the eddy turnover time need not be enforced, but is an output of the optimization process. Following this approach, two preferred spanwise wavelengths are found for the amplification of streamwise uniform streaks, for sufficiently high Reynolds numbers. The first ͑ z Ϸ 3h to 4h for plane channel flow͒ scales in outer units ͑i.e., with the channel half-width h͒, and the maximum amplifications increase with the Reynolds number. The second preferred spanwise length, which scales in inner ͑wall͒ units, is found at + Ϸ 100, in good accordance with the scales of the buffer layer cycles. Here milder energy amplification of initial conditions is found. These results strongly suggest that the mean lift-up effect is involved as an essential part of the mechanisms that permit the sustainance of turbulent coherent structures.

When streaks are of finite amplitude they may support the amplification of secondary perturbations ͓30,31͔ that lead to their breakdown, and to the refueling of streamwise vorticity necessary to sustain the streaks ͓32,33͔. This self-sustained process is christened in the emergence of nonlinear saddle solutions ͓34,35͔ and unstable periodic orbits ͓36͔. Many unstable nonlinear solutions have also been found for pipe flow ͓37,38͔, where a positive Lyapunov exponent is displayed in the turbulent regime ͓39͔. These "exact" solutions ͓40͔, which are believed to provide a structure to the subcritical transition process ͓41͔, typically have few unstable eigenvalues and associated modes, while they are attracting in all other directions in phase space. It has therefore been conjectured that transitional flow spends significant periods of time in the neighborhood of a few of these nonlinear solutions, and that therefore the lowest order statistics could be retrieved by averaging over the exact solutions. This rational framework has been boosted by the experimental observation of the velocity field bearing resemblance to exact nonlinear traveling wave solutions ͓42͔. Exciting phase space portraits of the transitional regime are becoming available for the plane Couette flow ͓43͔.

It is currently not clear how to extend the dynamical systems approach into the régime of developed turbulence. At large Reynolds numbers, exact solutions typically have high wave numbers, implying an unrealistic number of repetitions of the structure in space. Progress is being made with respect to this spatial dependence, with the discovery of localized counterparts to the infinitely extended solutions ͓44͔. Other periodic solutions have been continued to high Reynolds numbers ͓45͔. A different but serious problem is that, while capturing the typical features of turbulence near the wall, exact solutions do not exhibit the large dissipation of turbulence in the central region of the flow. It is this dissipation that leads to flattened mean profiles of the flow and strong deceleration of the core region for pipe flow. This slow core is clearly present in the observations of ͓42͔. While dissipation occurs on small scales in the central region, however, coherent streaks certainly are present. A myriad of solutions exhibiting streaks is now available for pipe flow ͓46͔, but it is not yet clear how to determine which best characterize turbulence. Toward this end, in the present study the eddy viscosity is invoked to help reveal on which length scales the most important structures occur, those associated with energy production through amplification, and who's presence may otherwise be masked by the dissipative eddies.

In addition to its potential relevance for understanding the structure of turbulent flows, the analysis of optimal energy amplifications is also relevant for flow control. In the laminar boundary layer, for instance, the forcing of nearly optimal streaks has been used to stabilize the base flow to Tollmien-Schlichting waves ͓47-49͔ and delay transition to turbulence ͓50͔. That large-scale coherent streaks can be effectively amplified in the turbulent boundary layer ͓51,52͔ has opened the way for their use in separation control in industrial applications ͓53͔. Furthermore, a few studies ͓54,55͔ have shown that by forcing large-scale streamwise vortices, viscous drag can be reduced in the plane channel flow. However, in these studies the cost of the control action has not been evaluated, and it is therefore unclear if a net power saving can be achieved. In this study, the advantage of forcing optimal structures is that the energy required for the forcing is minimized, reducing the cost of the control.

The scope of the present investigations is therefore to answer the following questions: what are the optimal energy amplifications and responses to forcing sustained by the mean turbulent flow in a pipe? Structures of which wavelengths are most amplified? And, how do these amplifications and their structure change with Reynolds number? Despite its ubiquity in practical applications, these results are not yet available for turbulent pipe flow. A linear model is presented in Sec. II. Linear results ͑Sec. III͒, based on the eddy-viscosity model, here and for the planar geometries, are yet to be compared with simulations of "real" turbulence, i.e., in full resolved simulation ͑Sec. IV͒. In other words, do the optimal streaks predicted by the linear eddy-viscosity model compare well to artificially forced streaks computed in direct numerical simulation? As simulations in the presence of turbulence are inherently nonlinear, what is the effect of finite amplitude forcing on these coherent streaks? Finally, we are interested in the influence of artificially forced finite amplitude streaks on the drag. Can the mean drag be reduced in the presence of large amplitude streaks? Is a net power saving, including the cost of the control, achievable, and if so, what is the best performance?

II. BACKGROUND

A. Turbulent base flow and eddy viscosity

Consider the incompressible flow of a viscous fluid of kinematic viscosity in a circular pipe of radius R. The bulk velocity U b ͑mean streamwise speed͒ is assumed constant. Lengths are nondimensionalized by R and velocities by U cl lam = 2U b , the center-line speed for laminar flow of the same bulk velocity ͓12,38,56,57͔. The usual Reynolds number based on the mass flux is defined as Re= 2U b R / . At sufficiently large Reynolds numbers the flow is turbulent and the turbulent mean flow profile preserves invariance in the spanwise direction , the streamwise direction x, and in time t. The mean flow velocity is written U = U͑y͒e x , where y = 1 r is the dimensionless distance from the boundary wall and r is the radial coordinate. In complement to this "outer" scaling, convenient for measuring large-scale properties, to describe near-wall structure it is commonplace to also define units via the friction velocity u 2 = ץ y U ͉ y=0 , based on the shear near the wall. The wall units of length and time are then / u and / u 2 , respectively. Variables nondimensionalized on this "inner" scaling bare the superscript + in the following, and the wall-Reynolds number is defined Re = u R / .

If, as a first approximation, the Boussinesq eddy viscosity is used to model the turbulent Reynolds-stresses, then the streamwise component of the Reynolds-averaged momentum conservation reads

1 Re ͩ 1 r + ץ r ͪ͑ T ץ r U͒ = ץ x P,

͑1͒

where the total effective viscosity is T ͑y͒ = 1 + E͑y͒ and E͑y͒ is the eddy-viscosity, normalized such that at the wall T ͑0͒ = 1, i.e., the kinematic value is attained. We denote by B = ץ- x P, the averaged streamwise pressure gradient necessary to maintain the prescribed mass flux.

For E͑y͒ we use the convenient expression originally suggested for pipe flow by Cess ͓58͔, later used for channel flows by Reynolds and Tiederman ͓59͔, then by many others ͓21,24,25͔:

E͑y͒ = 1 2 ͭ1+ 2 R ˆ2B 9 ͓2y -y 2 ͔ 2 ͑3 -4y + 2y 2 ͒ 2 ϫͫ1 -expͩ -yR ˆͱB Â+ ͪͬ 2 ͮ 1/2 - 1 2 . ͑2͒
The mean streamwise velocity U͑y͒ in equilibrium with this radially dependent eddy viscosity ͓Eq. ͑2͔͒ is easily inverted from the mean averaged momentum Eq. ͑1͒, and matches very well experimental observations. Here, R ˆ= Re/ 2, B ˆ= 2B. The fitting parameters A + = 27 and = 0.42 have been updated to improve the match with recent observations in ͓60͔.

B. Equations for the small coherent perturbations to the turbulent base flow

We now take the state ⌽ = ͑U , P͒, solution of Eq. ͑1͒, as a ͑steady͒ base flow and consider small coherent ͑correlated when ensemble-averaged͒ perturbations, = ͑u , p͒, to ⌽ in the presence of a coherent ͑correlated͒ force f and/or coherent initial conditions u 0 . Note that ⌽ is understood to be the turbulent mean state for f = 0 and u 0 = 0.

The equations governing the coherent perturbations require modeling of the Reynolds stresses to be closed. Reynolds and Hussain ͓29͔ discussed the eddy viscosity as a simple closure model, finding that it compared well for small-amplitude perturbations. Linearizing the closed equations, following ͓24-29͔, the equation governing smallamplitude coherent perturbations is

ץ t u + u r ץ r Ue x + Uץ x u = -١p + 1 Re ١ • ͓ T ͑yٌ͒͑u + ٌu T ͔͒ + f. ͑3͒
These perturbations also satisfy continuity, ١ • u = 0. The explicit components for the viscous term are given in Appendix A.

Exploiting the rotational and the streamwise translational invariance of the system, one may expand the perturbations in Fourier modes ͑r , , x ; t͒ = ͚ ␣m ˆ␣m ͑r ; t͒e i͑␣x+m͒ , where ␣ and m are the streamwise and azimuthal wave number respectively. Substitution into Eq. ͑3͒ gives ץ t ˆ␣m ͑r;t͒ = L ␣m ͑⌽͒ ˆ␣m + f ˆ␣m ͑r;t͒, ͑4͒

where L ␣m is the linear operator acting on the Fourier mode ˆ␣m ͑r͒. The problem may be considered for each mode independently, and the subscript ␣ , m is dropped in the following. The spanwise and streamwise wavelengths for each perturbation are given by = 2 / m and x = 2 / ␣ respectively in units of R.

C. Optimal growth and response to forcing

Three properties of the linearized system are now examined. They are the optimised growth of initial conditions, the optimised response to harmonic forcing, and the response to stochastic forcing. Standard definitions ͑see, e.g., ͓6,61͔͒ are briefly restated below for the present context.

Optimal growth in time over all possible nonzero initial conditions u ˆ0, for a particular ␣ , m mode, is given by the function

G͑␣,m;t͒ ϵ max u ˆ0 0 ʈu ˆ͑t͒ʈ 2 ʈu ˆ0ʈ 2 , ͑5͒
and the maximum achieved over all times is G max ͑␣ , m͒ ϵ max t G͑␣ , m ; t͒. Here f ˆ= 0, and the energy norm used is

ʈu ˆʈ2 ϵ ͑1 / V͒͐ V u • udV.
For harmonic forcing, f ˆ= f ˜ei⍀ f t , and accompanying response, u ˆ= u ˜ei⍀ f t , the optimal response is given by

R͑␣,m;⍀ f ͒ ϵ max f ˜ 0 ʈu ˜ʈ2 ʈf ˜ʈ2 . ͑6͒
The maximum over all frequencies

⍀ f is R max ͑␣ , m͒ = max ⍀ f R͑␣ , m ; ⍀ f ͒.
Finally, the response to a stochastic force is usually considered in the discretized form ͑here radially͒. Let f ˆ͑t͒ be the discretized stochastic force vector, and ͗ • ͘ denote ensemble averaging. Following analyses ͓9,22,61-64͔, the forcing is assumed to have Gaussian probability density with zero average ͑i.e., ͗f ˆ͘ = 0͒, isotropic ͑equal variance on all the three components͒ and delta-correlated in space and time, ͗f ˆ͑t͒f ˆH͑tЈ͒͘ = I␦͑t -tЈ͒, where the superscript H denotes the Hermitian transpose. The amplification of the stochastic forcing is then measured by the variance of the response:

V͑␣,m͒ ϵ tr͕C ϱ ͖, ͑7͒
where C ϱ = ͗u ˆu ˆH͘ is the covariance matrix for the response in the limit t → ϱ. This covariance matrix is computed by solving an algebraic Lyapunov equation ͑see Appendix B͒. The eigenfunction decomposition of the Hermitian matrix C ϱ is known as the Karhunen-Loève ͑KL͒ or proper orthogonal decomposition ͑POD͒. The ratio of its real and positive eigenvalues, n , to the variance V = ͚ n n represent the relative contribution to the variance from each orthogonal mode. The eigenfunction corresponding to the leading eigenvalue thus contributes most to the variance V. For details see ͓9,22,62-65͔.

D. Numerical methodology

The linear operator appearing in Eq. ͑3͒ must be discretized in order to compute the optimal growths and the responses described above. Some attention must be paid to difficulties introduced by the cylindrical geometry, and full details of the methods are given in Appendix B. Here we note the general approach and convergence of the results.

A Chebyshev-collocation method on up to N = 250 radial points was used to calculate up to 2N -3 eigenvectors and eigenvalues for the linear system ͓Eq. ͑3͔͒. The optimal growth code, which requires only these eigenvectors, eigenvalues and weights for calculating norms, is a well validated code used for several previous studies ͓25,27,28,66͔. Results for the optimal growth and the response to steady forcing were also verified by direct time stepping. At the largest Re= 10 6 the power spectral drop-off for the optimal modes was of eight orders of magnitude.

Responses from the linear model are compared with direct numerical simulations of the original Navier-Stokes equations. The pipe flow code described in ͓67͔ has been used for simulation of full turbulence subject to a forcing f. This time stepping code also uses a Fourier decomposition in and x, but finite differences on a nonequispaced radial ͑y͒ ninepoint stencil. Time steps are controlled using information from a second-order predictor-corrector method. For most simulations presented, the radial spacings are ␦r min + = 0.11, ␦r max + = 4.4 ͑wall units / u ͒. With dealiasing, ␦ + = 5.9, ␦ x + = 9.4 at the boundary, in a domain L = 15͑R͒, L + = 2700. Verification of these results is performed by calculation where the spacings above are respectively 0.08, 4.0, 3.9, 6.3 for L = 20. Simulations enforce a fixed mass flux, consistent with previous analyses ͓54,68͔.

III. LINEAR AMPLIFICATION OF COHERENT STRUCTURES

A. Wave number dependence

We first consider the dependence of the optimal amplifications on the streamwise and azimuthal wave numbers at the very large flux-Reynolds number of Re= 10 6 corresponding to Re = 19 200. Figure 1 shows the dependence of G max , R max and V on m for the selected ␣ = 0 , 1 , 10, 100, 1000. It is clear from all plots that streamwise elongated structures, where x Ͼ , are significantly amplified. The greatest growths are found for axially independent modes ͑␣ = 0͒, for which the curve provides an envelope over the results for nonzero ␣. The maxima on the ␣ = 0 curves are for m = 1 azimuthal symmetry. The inset to Fig. 1͑b͒ shows that steady forcing is optimal for axially independent harmonic forcing, ␣ = 0, while for streamwise nonuniform perturbations, a nonzero forcing frequency is optimal. For the stochastic forcing at the optimal value ␣ = 0 , m = 1, it is found that the leading Karhunen-Loève mode dominates, representing 93% of the contribution to the total variance, as reported in the inset to Fig. 1͑c͒.

A secondary peak occurs in G max at a larger m corresponding to an azimuthal wavelength of + = 92 ͑wall units / u ͒. Contrary to the optimal growth of initial conditions,

n α = 0, m = 1 λ + θ (a) (b) (c)
FIG. 1. ͑Color online͒ Optimal energy growth of initial conditions G max ͑a͒, optimal response to harmonic forcing R max ͑b͒ and variance of the response to stochastic forcing V ͑c͒ for Re = 19200. Dependence on azimuthal wave number m, with corresponding wavelength in wall units + ͑top axis͒, for selected values of the dimensionless streamwise wave number ␣. Reference slopes m -2 and m -1 ͑dotted lines͒ are shown in subfigures ͑b͒ and ͑c͒ respectively. The dependence of the optimal harmonic forcing amplification R͑␣ , m = 1 ; ⍀ f ͒ on the forcing frequency ⍀ f is reported in the insert of subfigure ͑b͒ for ␣ = 0 and ␣ = 1; R max = max ⍀ f R. The percentage contribution to the variance, 100 j / V͑0 , 1͒, of the response to stochastic forcing contained in the first 20 Karhunen-Loève modes is reported in the insert of subfigure ͑c͒.

no secondary peak is found in the optimal response to forcing. A closer examination of Figs. 1͑b͒ and 1͑c͒ reveals that the R max and V curves corresponding to ␣ = 0 scale approximatively like m -2 and m -1 , respectively, for intermediate values of m. A noticeable change of slope is found for values of m corresponding to + Ϸ 100. In plane channel flow, this power-law dependence has been explicitly related to the behavior of geometrically similar optimal structures in the loglayer ͓26͔. In the case of the response to forcing, primary and secondary peaks occur in the deviations from this power-law behavior ͓26͔. This can be seen by considering the premultiplied responses m 2 R max and mV. A secondary peak at + = 73 is found in the premultiplied response curves for ␣ = 0, shown in Figs. 2͑b͒ and2͑c͒.

B. Reynolds number dependence

The spanwise optimal wave numbers, corresponding to the primary and to the secondary peak, do not change with Re when, respectively, scaled with outer and inner units. This can be seen in Fig. 2, where G max and the premultiplied m 2 R max , mV computed for Re= 10 4 , 10 5 , and 10 6 ͑corresponding to Re = 317, 2380, and 19 200, respectively͒ for ␣ = 0 are reported. While the premultiplied response to steady forcing is useful for showing deviations from the R max ϳ m -2 trend, the premultiplied stochastically forced case is arguably more representative of turbulence. When m is presented on a log scale, equal area under the curve mV implies equal contribution to the total variance. Here, where the spanwise wave number is integer for pipe flow, all the primary peaks are reached for m = 1 except for the premultiplied stochastic response at lower Re ; for Re = 10 4 , m = 2 is preferred and the selection of m is very weak in mV for Re = 10 4 and lower. The secondary peak is always found at + Ϸ 70-90 for all the premultiplied responses. The primary and the secondary peaks, however, are not yet separated at Re = 10 4 , while they are for Re = 10 5 and higher, in good accordance with the turbulent plane channel ͓24-26͔ the boundary layer ͓27͔.

The maximum responses, associated with the primary peak ͑␣ = 0, m = 1͒ increase with the Reynolds number, while those associated with the secondary peak remain almost constant. The structure of the leading modes corresponding to

+ m + V + 10 -5 (a) (b) (c) (d) (e) (f)
FIG. 2. ͑Color online͒ Dependence on the azimuthal wave number m and the corresponding azimuthal wavelength ͑top axis͒ of: the optimal energy amplification of initial condition G max ͓͑a͒ and ͑d͔͒, the premultiplied optimal response to harmonic forcing m 2 R max ͓͑b͒ and ͑e͔͒ and of the premultiplied variance of the response stochastic forcing mV ͓͑b͒ and ͑e͔͒. Streamwise uniform ͑␣ = 0͒ perturbations are considered for the selected Reynolds numbers Re= 10 4 , 10 5 and 10 6 corresponding to Re = 317, 2830, and 19200. The responses, wave numbers and wavelengths are scaled in outer units on the upper row ͓͑a͒, ͑b͒, and ͑c͔͒ and in inner units in the lower row ͓͑d͒, ͑e͒, and ͑f͔͒. the primary peak for optimal growth, for the response to steady and stochastic forcing, appear to be almost identical and remain the same over the entire range of Reynolds numbers. Plotted in Fig. 3 are the leading modes for m = 1 , 2 , 4 ͑data from the optimal harmonic forcing case at Re= 10 6 ͒. The optimal input consists of the expected large-scale streamwise vortices and the optimal output of large-scale streamwise streaks. Note that in this geometry, only modes of m = 1 can exhibit flow across the axis. The structures corresponding to the secondary peak, which scale in inner units, are almost indistinguishable from those found in the plane channel and the boundary layer ͓25-27͔ and are therefore not reported here.

IV. SIMULATION OF FORCED FINITE AMPLITUDE COHERENT STREAKS

A. Response to finite amplitude optimal steady forcing

The optimal response predictions, described in the previous section, were based on two main simplifying assumptions: ͑a͒ the effect of the uncorrelated motions was modeled with an isotropic eddy-viscosity and ͑b͒ coherent perturbations were assumed small. We now test the level of approxi-mation implied in these assumptions by computing in full DNS the effective responses of the predicted optimal perturbations. In particular, we consider the response to the optimal steady large-scale forcing, the mode which is the most amplified, and potentially the most relevant for passive flow control applications. The forcing term f is taken to be proportional to the mode f ˜that optimizes R max , obtained from the above linear analysis, for ␣ = 0 and a given m. As the considered optimal forcing is steady, time averages can be easily calculated and the statistical convergence readily checked. The magnitude of the steady forcing ʈfʈ is measured by the rms value, as before. The response of the flow in simulations is given by u ¯͑r , ͒ = ͗u͘ x,t -U ¯͑r͒x ˆ, where U ¯͑r͒ = ͗u͘ ,x,t is the mean flow in the absence of forcing. ͗ • ͘ s indicates averaging over the subscripted variables. Averages are taken over at least 2800R / U cl lam time units. The DNSs are performed at Re= 5300, for which Re = 180.

From Fig. 4 it is seen that the amplitude of the responses ʈu ¯ʈ depends linearly on the forcing for small forcing amplitudes ʈfʈ. Small forcing can lead to large responses, in agreement with the linear eddy-viscosity model and, in the linear régime, the leading mode is m = 1, as expected. The response predicted by the linear eddy-viscosity model is R max = 276 2 for m = 1 at this Re, while in the simulations, ʈu ¯ʈ / ʈfʈ is a factor 2 smaller than expected in the linear range. As the time-averaged mean flow for the unforced case, obtained from DNS, deviates very little from that given by Eqs. ͑2͒ and ͑1͒, this difference is most likely attributed to the isotropic eddy-viscosity assumption. Similar considerations apply to the observation that the m = 4 mode is slightly more responsive than the m = 2 mode in the linear regime.

For larger forcing amplitude, the amplitude of response saturates. For the m = 1 mode, this saturation begins when the amplitude A S of the associated coherent streaks is near to 19% of the center-line speed, where A S = 1 2 ͑max u ¯x -min u ¯x͒ ͓69͔. The responses to higher order optimal modes ͑m = 2 and m = 4͒ saturate at larger forcing amplitudes ʈfʈ, and therefore contain more energy ʈu ¯ʈ2 than the m = 1 mode, for sufficiently large forcing amplitudes.

In Fig. 5, we report the finite amplitude coherent vortices and streaks induced by the steady forcing for, respectively, low and moderate forcing amplitudes for the three azimuthal modes m = 1, m = 2, and m = 4. From these figures, and by comparing them with the linear optimals reported in Fig. 3, it is seen how, for increasing forcing amplitude, the low-speed streaks are more and more spatially concentrated, while high-speed streaks have a tendency to widen.

B. Influence on the drag

We now examine the influence of the forced finite amplitude streaks on the drag, always at Re= 5300 ͑corresponding to Re = 180͒. Figure 6 shows time series for the instantaneous C f = 2͑u / U b ͒ 2 in full DNS for three levels of optimal steady forcing of the m = 1, ␣ = 0 mode. In wall units, the spanwise wavelength of the forced streaks for m = 1 is + Ϸ 1130, which extends by a factor greater than 2 the spacings considered in earlier investigations of the plane channel flow ͓54,55͔. The horizontal line in Fig. 6 represents the baseline ͑time averaged͒ C f in the absence of forcing. From this figure we see that indeed average drag reductions can be obtained also in the turbulent pipe flow. We have verified the accuracy of this result, in particular that no loss of the observed drag occurs ͑solid line͒, by performing an additional calculation at higher resolution ͑more than doubling the number of grid points, see Sec. II D͒.

In Fig. 7, the relative change for time averaged C f is reported versus the forcing amplitude for the m = 1, m = 2, and m = 4 modes. Drag reduction is obtainable over a range of forcing amplitudes covering approximately one order of magnitude. Selecting the best forcing amplitude, the maximum drag reduction achieved is 12.8% when forcing the m = 1 mode. It is reduced for higher modes, e.g., 10.4% for m = 2 ͑ + Ϸ 565͒ and 4% for the mode m = 4 ͑ + Ϸ 280͒. Forcing streaks at the near-wall peak mean spacing ͑ + = 100, corresponding to m Ϸ 12͒ results in increased drag.

For this control strategy to be interesting for applications, the net power saved must be considered. The power con-sumption ͑per unit length͒ spent to maintain the desired constant volume flux of the turbulent flow is P = 8 C f in units U cl lam3 R, and the control power spent ͑per unit length͒ to maintain the vortices is W = 1 L ͐u • fdV. The drag relative to the unforced case is then given by C f / C f0 = P / P 0 , and the total relative net power consumption by ͑P + W͒ / P 0 . The percentage net power saving is thus 100͓1 -͑P + W͒ / P 0 ͔. From Fig. 7, where the relative drag and net power consumption are shown versus the forcing power, it is seen that the best net power saving corresponds to 11.1% and is obtained for the m = 1 mode. For this case, W is only 1.4% of P 0 . The response to forcing the large-scale mode is large, and as a result little energy is required to generate the prominent large-scale slow streak which affects the drag at the walls. The forcing of higher optimal modes results in smaller net power savings ͑e.g., at best 8.3% for m = 2 and 3.3% for m = 4͒.

The drag reducing mechanism appears to be consistent with that observed in plane channel flow ͓54,55͔. Figure 8͑a͒ shows a snapshot of streamwise vorticity field in the absence of forced streaks. Intense small-scale quasistreamwise vortices with + Ϸ 100 appear irregularly but in the whole azimuthal domain. These vortices, largely responsible for the turbulent skin friction ͓70͔, are weakened when large-scale streaks are artificially forced, and are almost suppressed from the high-speed streak region as seen in Figs. 8͑b͒ and8͑c͒. In the absence of this mechanism the drag would only increase due to distortion of the mean flow by the forcing. This drag, however, only appears to be large for moderate to large amplitudes of the forcing, leading to the observed drag reduction for intermediate forcing.

V. SUMMARY AND CONCLUSIONS

The first part of this investigation was dedicated to the computation of the optimal linear energy amplifications of coherent structures supported by the turbulent pipe mean flow. The problem formulation and the tools used in the analysis follow those of investigations of the turbulent plane channel flow ͓24-26͔, Couette flow ͓28͔, and boundary layer flow ͓27͔. In accordance with these investigations it is found that, despite the turbulent mean flow is linearly stable, it can support the energy amplification of initial conditions and of harmonic and stochastic forcing. These energy amplifications have the expected characteristics of the lift-up effect, where vortices amplify streaks and the most amplified perturbations are streamwise uniform ͑␣ = 0͒. For these streamwise uniform structures, at sufficiently large Re, the maximum energy amplifications are obtained for the m = 1 azimuthal mode. Furthermore, this large amplification increases with the Reynolds number. As observed in the Couette and plane channel flow ͓26,28͔, for the three different calculations, the initial value problem, response to harmonic forcing and most energetic Karhunen-Loève mode for stochastic forcing, the optimal inputs and outputs almost coincide.

A secondary peak is found in the G max ͑m , ␣ = 0͒ curves for the azimuthal mode corresponding to an azimuthal wave number + = 92, but only a change of slope in the R max ͑m , ␣ = 0͒, V͑m , ␣ = 0͒ curves is noticeable at this + . The responses R max and V of optimal harmonic and stochastic forcing, for ␣ = 0 decrease with the azimuthal wave number, with an approximate scaling of m -2 and m -1 respectively for intermediate m. For the plane channel flow, this inverse power dependence on the spanwise wave number has been related to the amplification of geometrically similar structures in the log-layer ͓26͔. To highlight the departure of the optimal responses from the power-law trend, the premultiplied responses m 2 R max and mV have been considered. A secondary peak is then found at + = 73. This peak wavelength and the corresponding responses do not change with Reynolds number, provided that they are scaled in inner units, and that the Reynolds number is sufficiently large, so as to allow for a well defined inner-outer scale separation. The corresponding structures consist of streamwise uniform buffer-layer streaks and vortices.

The linear formulations we have used allowed a relatively straightforward computation of the linear optimal energy amplifications. The advantage of this model, based on the turbulent eddy viscosity, is therefore its simplicity, while quantitatively accurate predictions could be made with more sophisticated models. It is however important to verify that the predictions of the model, such as the order of magnitude of the growth of structures at large scale, or the fact that they more amplified than small-scale structures, are at least qualitatively verified. We have therefore reported the results of direct numerical simulations, where coherent structures have been artificially forced at the moderate Re = 180. The optimal steady forcing of large-scale structures ͑m = 1, m = 2, and m = 4͒ has been considered as these structures are the most responsive, and therefore are not only easier to detect, but are potentially most interesting for control purposes. The main predictions of the linear model are qualitatively observed in the DNS: the low-m optimal perturbations show a very large response in the amplitude of induced streaks. For sufficiently small forcing amplitudes, the response to forcing is linear and the azimuthal symmetries of the forced modes are well preserved in the response. In this linear regime, the m = 1 mode is the most responsive. The agreement between the amplifications predicted by the linear model and the one observed in the DNS is only qualitative. The quantitative differences can probably be attributed to the very eddyviscosity assumption and could probably be reduced by using a more realistic modeling of the Reynolds stresses. As coherent structures are strongly elongated in the streamwise direction, the isotropic eddy viscosity is gives only a first order approximation. At finite amplitudes of the coherent motions, the neglected interactions with the mean flow and the random turbulent field clearly become important ͓71͔. For larger forcing amplitudes, a saturation in the amplitude of the forced response is observed. The m = 1 mode is the first to saturate, at a considerable amplitude of almost 20% of the centerline velocity. For increasing amplitudes, the symmetry between the high and low-speed streaks is progressively lost. Low speed streaks concentrate in smaller spanwise extensions, with slightly larger amplitude, while highspeed streaks widen and slightly decrease in amplitude. This asymmetry recalls the observed "real" unforced turbulent flows ͓72͔.

The exact relation existing between optimally amplified streaks and the coherent structures populating the "natural" ͑unforced͒ turbulent pipe flow requires consideration. In the plane channel and Couette flows, the similarity of the spanwise scales selected by the linear mechanism and those observed in the actual turbulent flow is probably explained by the requirement that significant amplification is necessary for energetic turbulent processes to be self-sustained. What is currently not clear, because of the lack of experimental or numerical data, is if m = 1 is the dominating coherent mode at sufficiently large Re. Most of the investigations of large-and very-large scale coherent motions in the pipe flow ͓19,20,73,74͔, have rather focused on the streamwise scales of these motions. If these motions arise through a selfsustaining process, then the streamwise scales of the motion would be selected by processes related to the breakdown of the streaks. This process is missing in the present analysis. What can be said with the model, however, is which azimuthal modes are most amplified for a streamwise scale selected by other stochastic processes. For instance, if one considers streamwise wavelengths ranging from x = 16͑R͒ ͑␣ Ϸ 0.4͒ to x = 6 ͑␣ Ϸ 1͒, then the most amplified structures vary from m = 2 to m = 5 modes at sufficiently large Re in the stochastic response. However, at even moderately low Reynolds number ͑e.g., Re= 10 4 corresponding to Re = 317͒, the azimuthal selection is weak and modes ranging from m = 1 to m = 10 all exhibit a comparable variance ͓see Fig. 2͑c͔͒ in the integrated response. These predictions are compatible with recent results from direct numerical simulations at Re = 150 ͓75͔ where the most energetic POD modes are found for the streamwise wavelength x = 20͑R͒ ͑the computational box size͒ and for azimuthal modes ranging from m = 2 to m = 7. Additional simulations or experiments would be highly desirable to examine the azimuthal structure of the most energetic modes at larger Re.

One of the most relevant applications of the optimal perturbations of the turbulent mean flows probably resides in flow control. The direct numerical simulations confirm that energy can be greatly amplified in a turbulent flow. An objective of both theoretical and practical relevance is the reduction of the power required to drive the constant mass flux through the pipe. Using optimal forcing, less than 2% of the power required to pump the ͑unforced͒ flow, is required to induce modifications in the mean flow of the order of 20% of the centerline mean velocity. Previous numerical ͓54͔ and experimental ͓55͔ studies have shown that it is possible to reduce the mean turbulent drag in the plane channel flow by forcing large-scale vortices with spanwise wavelength of Ϸ400-500 wall units. Here we have extended this result to pipe flow, where + Ϸ 1200 wall units for m = 1. Mean drag reductions of Ϸ13% are achieved by forcing this m = 1 mode, while the drag reduction is reduced for reduced spanwise spacing. In ͓54͔ the observed drag reduction was attributed to the mean weakening of the circulation of quasistreamwise coherent vortices. This effect on the quasistreamwise vortices is found in the present turbulent pipe flow. Our visualizations show that the quasistreamwise vortices of the "natural" buffer layer cycle are packed in the low-speed streak region, while they are essentially removed from the highspeed streak region.

At present it is not clear why the drag reduction reduces as the spanwise spacing is reduced. It is possible that for m = 2 ͑ + Ϸ 600͒ the reduced scale separation from that of the naturally occurring streaks ͑ + Ϸ 100͒ causes the dragreduction effect to be reduced. At spanwise separations greater than those possible at this Re , where clear scale separation occurs, a plateau in the drag reduction is possible for separations greater than + Ϸ 1200. However, if the greatest reduction continues to be observed for the largest possible spacing, here m = 1, it may be conjectured that, rather than scaling in inner units, the mode that most reduces drag has an upper limit that depends on the geometry. Additional simulations or new experiments at much larger Re could confirm or reject these conjectures.

We have also extended earlier investigations by computing the net power saving obtained by forcing the optimal coherent streaks. As these streaks are optimal, the power required to force them is already almost minimized. Even at the considered low Reynolds number, the controlling streaks are forced with less than 2% of the power necessary to maintain the base flow, and a net power saving of Ϸ11% is achieved with the m = 1 mode. For useful implementation of the method proposed in this study, the induction of rolls need not be efficient. For example, efficiency of 50%, still using only 4% of the driving energy, might be acceptable. Note also that, as the amplification of the forcing is expected to increase with increasing Re, the percentage of total power used to force the streaks is predicted to decrease for increasing Re. In practice, induction of large-scale rolls is possible via passive actuators in laminar ͓49,50,76͔ and turbulent flows ͓51-53͔. Although neither would reproduce the precise details of forcing used here, linear optimals are known to be robust to perturbations of the system, and responses close to those seen here can be expected. An experimental verification of the proposed control strategy would of course be very welcome.

FIG. 3 .

 3 FIG.3. ͑Color online͒ Linear optimal responses to steady forcing for ␣ = 0 and m = 1 , 2 , 4; Re = 19200. Vectors: cross-stream components of the input vector field. ͑Colored͒ contour levels: streamwise component of the output field ͑streaks͒. White fast-streaks, red/dark slow streaks.

  FIG. 4. ͑Color online͒ Response to forcing in full simulationshowing the linear régime and nonlinear saturation. ʈfʈ is the rms value. The isolated line is of slope 1 and squares correspond to the plots of Fig.5.

  FIG. 5. ͑Color online͒ Nonlinear responses to finite amplitude forcing for ␣ = 0 and m = 1 , 2 , 4 from direct numerical simulation at Re = 180. Small forcing amplitude ͑upper row͒ and moderate forcing ͑lower row͒, corresponding to the marked squares in Fig. 5. ͑a͒ m = 1, ʈfʈ = 8.010 -5 ; ͑b͒ m = 2, ʈfʈ = 1.410 -5 ; ͑c͒ m = 4, ʈfʈ = 8.010 -5 ; ͑d͒ m = 1, ʈfʈ = 8.010 -4 ; ͑e͒ m = 2, ʈfʈ = 8.010 -4 ; ͑f͒ m = 4, ʈfʈ = 2.510 -4 . Vectors: cross-stream components of the forced velocity field. ͑Colored͒ contour levels: streamwise component of the output field ͑streaks͒. White fast-streaks, red/dark slow-streaks.

  FIG.6. ͑Color online͒ Relative drag for several levels of steady forcing of the m = 1 mode from direct numerical simulation. At intermediate levels of forcing the drag is reduced, falling below the average of the unforced case ͑horizontal straight line͒. The forcing can reduce the drag. The solid line is a higher resolution check at ʈfʈ = 8.010 -4 .

  FIG.8. ͑Color online͒ Snapshots of the cross-stream distribution of the streamwise vorticity x field at a selected streamwise station from direct numerical simulation. ͑a͒ Streamwise vorticity in the absence of forcing. ͑b,c͒ Weakened and localized streamwise vorticity in the presence of the m = 1 , 2 steady forcing, corresponding to optimal drag reduction. Blue/dark intense vorticity, white no vorticity. Same color scales are used on all the subplots. Contour lines denote isolevels of the streamwise averaged velocity ͑the forced streaks͒.

ACKNOWLEDGMENTS

We gratefully acknowledge financial support from the European Community's Seventh Framework Programme ͑FP7/ 2007-2013͒ under Grant agreement No. PIEF-GA-2008-219-233.

APPENDIX A: RADIALLY-DEPENDENT VISCOSITY

For the radially dependent effective viscosity, T ͑r͒, the viscous terms in Eq. ͑3͒ are given by

͑A1͒

where

APPENDIX B: NUMERICAL CONSIDERATIONS

Cumbersome differential operations in cylindrical coordinates, particularly given T = T ͑r͒, mean that finding the Orr-Somerfield form for Eq. ͑B2͒ is not straight forward in this geometry. Instead, for this study the continuity equation together with Eq. ͑3͒ have been used directly in primitive variable form, giving the generalized eigenvalue problem

where ͑r ; t͒ = ˆ͑r͒e t . Note that Y is singular, containing three zero boundary conditions and a zero right hand side for the continuity condition. Simplifications can be made by observing that X is invertible, and that, as the have real part strictly less than zero, Eq. ͑B1͒ can be written in the standard form, ͑-1 / ͒ ˆ= ͑X -1 Y͒ ˆ. Standard LAPACK routines re- turn the leading eigenvalues along with a predicable number of infinite eigenvalues that are easily filtered. Excellent numerical stability is observed.

The continuity equation together with the linearized Eq. ͑3͒ and no-slip boundary conditions can be written in the form u = Vx,

where x is an alternative representation for the state space .

Conversion between the representations is performed by the action of V and W. This form is again familiar through the Orr-Somerfield formulation of the governing equations.

Here, however, the eigenfunction decomposition is employed. Matrix V has columns of eigenvectors u ˆn and ⌳ = diag͑ n ͒, with the n arranged in descending order. These matrices, along with the diagonal weight matrix M s.t. ʈuʈ 2 ϵ u H Mu, are sufficient to calculate optimally growing initial conditions. Consider the norm for x given by ʈuʈ 2 = x H M ˆx ϵ ʈxʈ 2 = ͑x , x͒, where M ˆ= V H MV. The optimal initial condition x 0 for growth at time t satisfies G͑t͒ = ͑Ax 0 , Ax 0 ͒ / ͑x 0 , x 0 ͒ = ͑A + Ax 0 , x 0 ͒ / ͑x 0 , x 0 ͒, where the transfer function for x is A = e t⌳ and the adjoint can be shown in a few operations to be A + = M ˆ-1 A H M ˆ. The growth and initial condition can be calculated by power iteration on the matrix A + A. This matrix is straight-forward to evaluate when ⌳ is diagonal.

Calculation of the optimal harmonic forcing requires only minor alterations to the above. In this case one uses A = ͑i⍀ f I -⌳͒ -1 , which is again a diagonal matrix for ⌳ diagonal.

For the stochastic response, it can be shown ͓65͔ that C ϱ = VG ϱ V H , where G ϱ solves the Lyapunov equation

For W one may take the Moore-Penrose pseudoinverse of V. ͑Singularity of Y leads to fewer columns of eigenvectors than rows in V.͒ MATLAB routines are then available to solve Eq. ͑B3͒.