Rabah Bouzidi

Sandrine Buytet

Anh Le Van

A numerical and experimental study of the quasi-static deployment of membrane tubes

The paper is concerned with the deployment of inflatable structures in quasi static conditions. The prob lem is formulated as a minimization problem under constraints with emphasis on the self contact aspect and it is solved by a descent algorithm of the first order which allows to overcome the singularity or ill conditioning difficulty specific to membrane structures. Numerical computations are carried out on a simply folded membrane tube and compared with experimental results in order to prove the ability of the proposed numerical approach to correctly predict the deployment process.

Introduction

Inflatable structures are membrane like structures which acquire their stiffness and bearing capacities when inflated by air or another gas. They possess highly interesting mechanical proper ties which account for their increasing development in space industry: compared to standard structures, they are light, easy to be folded, rapidly deployable, their storage volumes are small, their manufacturing is less expensive, and moreover they have better thermal properties. The successful inflatable antenna experiment (IAE) carried out by [START_REF] Freeland | Significance of the inflatable antenna experiment technology[END_REF] showed that inflat able structures can be used in space applications as an alternative to conventional electromechanical deployable systems. Many types of inflatable structures have then been developed especially in building large scale systems such as solar arrays, solar sails, ra dar antennas, thermal or light shields, and a lot of studies have been devoted to them. Yet not everything is easy with the deploy ment of inflatable structures. The main difficulty is to control the change of geometry all along the deployment stage and further experimental or numerical investigations are still necessary in order to better understand the deployment mechanism.

As far as the numerical computations are concerned, one has to take into account different types of nonlinearities due to the large deformations and finite strains, the material nonlinearities, the inflating pressure which is a follower loading, the contact of the membranes with external obstacles and particularly the self contact. Bifurcation may also be involved as wrinkles and folds can appear in the membranes. Such nonlinear problems are usually dealt with by means of the standard finite element method based on the discretization of principle of virtual work, where the equa tion of motion is discretized and solved at each increment by an iterative Newton like scheme. The specific difficulty encountered is due to the fact that the in plane stiffness is much higher than the flexural stiffness which is quite negligible, and, as a conse quence, the stiffness matrix is almost singular as long as the inflatable structure is not fully deployed and taut. This singularity can be avoided by considering the deployment in the dynamic context, as done by [START_REF] Salama | Simulation of the deployment dynamics of inflatable structures[END_REF], [START_REF] Clem | A pressurized deployment model for inflatable space structures[END_REF] and [START_REF] Wang | Deployment simulation of ultra-lightweight inflatable structures[END_REF] using an explicit scheme. In the first computational steps when the stiffness matrix is singular, the inertial terms allow the global tangent matrix to be regularized. In general, the deployment scenario depends on the rate of the inflating pressure, as shown in [START_REF] Wang | Deployment simulation of ultra-lightweight inflatable structures[END_REF].

On the other hand, the dynamical solution may be difficult especially when the deployment takes place in so low speed that the inertial forces are not significant. In this case, one had better abandon the dynamic point of view and solve the deployment problem in the quasi static framework. Few works have been devoted to the self contact of the inflatable membranes in large deformations. The aim of this paper is to propose a numerical approach to solve the quasi static deployment problem, which cope with the above mentioned singularity difficulty in a simple and efficient way. The self contact is assumed to be frictionless so that all forces are conservative and the solution is obtained by minimizing a convenient total potential energy. A descent algo rithm of first order will be chosen to achieve the minimization, as done in [START_REF] Bouzidi | Numerical solution of hyperelastic membranes by energy minimization[END_REF] and [START_REF] Gil | Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains[END_REF]. Although first order methods have only a linear rate of convergence, they have a great advantage over the Newton iterative scheme of requiring the gradient of the energy only not the stiffness matrix thus allow ing one to compute all the equilibrium configurations, even when the inflatable structure is not fully stretched. The proposed ap proach is particularly interesting in the case of weak inflating pres sures, as those used in space applications.

The paper is outlined as follows. First, the basic equations of the contact problem will be set out and then it will be shown how to recast the resulting strong formulation into a minimization prob lem under constraints. Next, the algorithm of the projected gradi ent method used to solve the minimization problem and the finite implementation will be described. Numerical computations of the deployment will be carried out on a simply folded membrane tube in Section 6. Experimental results obtained on the same tube will be presented in Section 7, which prove the ability of the proposed numerical approach to correctly predict the deployment process.

Strong formulation of the self-contact problem

The membrane structure

Let us consider a membrane structure with no bending stiffness, which undergoes large deformations so that it may come into self contact over a part of its surface and during some time interval. We shall describe the basic equations of the problem using a Lagrang ian description. The reference position of the membrane is defined by the reference middle surface S 0 and the reference thickness H. The current position S of the middle surface at any time t in the time interval ½O; T of interest is given by the deformation /: 8t 2 ½O; T; /ð:; tÞ : S 0 3 X#x /ðX; tÞ ð 1Þ where X and x are the reference and current positions of any parti cle on the middle surface. The positions X and x are parameterized by the curvilinear coordinates n ðn 1 ; n 2 Þ via the following chain (Fig. 1):

A region of R 2 ! h S 0 ! / S n ðn 1 ; n 2 Þ # X hðnÞ # x /ðX; tÞ jðn; tÞ ðj / hÞ ð2Þ
Assuming that all the mappings introduced are smooth enough, the covariant local base on the reference middle surface is defined as

A a X ;a @h @n a A 3 N A 1 Â A 2 kA 1 Â A 2 k ð3Þ
where Greek indices are equal to 1 or 2; A 3 N is the unit normal vector at point X 2 S 0 . Similarly, the covariant local base on the current middle surface is defined as a a x ;a @j @n a n a 3

a 1 Â a 2 ka 1 Â a 2 k ð4Þ
The covariant components of the Green strain tensor E on the middle surface are

E ab 1 2 ða ab A ab Þ ð 5Þ
where A ab A a :A b and a ab a a a b are the covariant components of the metric tensor in the reference and current configurations, respectively.

The membrane is assumed to be made of a hyperelastic material defined by a strain energy density wðEÞ per unit volume. The sec ond Piola Kirchhoff stress tensor R is related to strain tensor E via the constitutive relationship

R

@wðEÞ @E ð6Þ

The plane stress assumption in the membranes is satisfied by zeroing the normal stress, R 33 0, in the energy density evalua tion. The reference surface S 0 is partitioned into three non overlap ping and pairwise disjoint parts denoted S 0U ; S 0T and S 0c , where S 0U and S 0T are the parts where the displacement U x X and the tractions are prescribed, respectively, and S 0c is the part where self contact potentially takes place during the time interval ½O; T. where U; T are the prescribed displacement and nominal traction, respectively; P is the first Piola Kirchhoff (nominal) stress tensor.

The boundary conditions on the contact surface are deferred to Section 2.4 where the contact law is specified.

Contact kinematics

Contact, or more precisely self contact, is one of the chief agents of the deployment and has to be taken into account in the analysis.

It is assumed that one portion S ð1Þ 0c of the reference membrane sur face can potentially come into contact with another portion S ð2Þ 0c (so that the above mentioned contact surface S 0c is S

0c [ S ð1Þ
The mappings h ðiÞ and j ðiÞ ð:; tÞ /ð:; tÞ h ðiÞ are specific to the contact surfaces, they are not necessarily the same as the membrane parametrization (2). The target surface S ð2Þ c should be conveniently parametrized so that the normal at point x ð2Þ to S ð2Þ c n ð2Þ @j ð2Þ @n ð2Þ1 Â @j ð2Þ @n ð2Þ2 @j ð2Þ

@n ð2Þ1 Â @j ð2Þ

@n ð2Þ2 ð9Þ

is directed toward the contactor S ð1Þ c . Now let x be a particular point on contactor surface S ð1Þ c , one defines the opposite point y 2 S ð2Þ c as the closest point to x in the Euclidean sense (Fig. 3):

y arg min x ð2Þ 2S ð2Þ c kx x ð2Þ k ð 10Þ
In writing (10), we assume that point y is unique without dis cussing the case of non uniqueness, which is not of major conse quence in numerical computations. The interested readers may consult comprehensive comments on this issue in the literature, e.g., [START_REF] Heegaard | An augmented Lagrangian method for discrete large-slip contact problems[END_REF], Wriggers (2002, p. 48), Laursen (2002, p. 115), and references quoted therein. If the target surface S ð2Þ c is smooth at y, then y is the projection of x on S ð2Þ c . As may be seen, the notation x ðiÞ ; i 2 f1; 2g, in (8) designates any point on where m is the normal at point y. The sign convention chosen in (11) entails that the proximity g is positive when there is interpenetra tion and negative when there is a gap between the two bodies. Thus, symbol g is the opposite of what is referred to as the gap in the literature. In order to check that g is an objective quantity, let us consider a change of reference frames defined by

x à Q ðtÞ Á x þ cðtÞ ð 12Þ
where vector cðtÞ represents the relative rigid body translation between the original reference frame and the new one (denoted by superscript ⁄ ) and Q ðtÞ is an orthogonal tensor representing the relative rotation of the reference frames. From relations

x à y à Q ðtÞ ðx yÞ and m à Q ðtÞ m, one derives

g à m à Á ðx à y Ã Þ Q ðtÞm ð Þ Á Q ðtÞðx yÞ ð Þ
m Á Q T ðtÞQ ðtÞ Á ðx yÞ g ð13Þ using equality Q T ðtÞQ ðtÞ I (identity tensor). Relation (13) shows that the proximity g is objective.

Let us now convert all the spatial quantities defined above into material ones defined on the reference configuration, by means of the deformation / in (1). Considering the point X 2 S ð1Þ 0c related to the given point x by x /ðX; tÞ, and the point Y 2 S ð2Þ 0c related to point y by y /ðY; tÞ, the inverse image of X (resp. Y) under the parametrization h ð1Þ (resp. h ð2Þ ) will be denoted n (resp. g). Points n; X and x are linked through the following chain:

A region of R 2 ! S ð1Þ
For notational conveniences, the dependencies on / and j ð2Þ will be omitted in the material kinematics variables and only the dependency on ðX; tÞ will be shown. Thus, for instance, one will write g gðX; tÞ and g gðX; tÞ.

Contact tractions

We assume frictionless contact during the deployment process and resolve the contact Cauchy traction vector tðx; tÞ at a typical point x 2 S ð1Þ c in terms of the normal vector m at the projection point

y 2 S ð2Þ c as follows t t N m ð17Þ
The normal component t N of traction vector t in the direction of normal m is the contact pressure, positive when the contact occurs.

Consider the change of reference frames defined by ( 12) again. Since the Cauchy traction vector t itself is objective, t à Q ðtÞ t, one can recast the contact pressure in the new frame as

t à N t à Á m à Q ðtÞ Á t ð Þ : Q ðtÞ Á m ð Þ t Á Q T ðtÞQ ðtÞ Á m t Á m t N ð18Þ
which means that the contact pressure t N is an objective variable.

Let TðX; tÞ P N be the nominal Piola Kirchhoff traction

vector at point X 2 S ð1Þ
0c (N is the normal vector at X to S ð1Þ

Contact law

The contact law describing the interfacial response at current time t is naturally expressed with spatial variables. Here, it is recast o into the mixed form involving the spatial proximity and the nominal contact tractions, which is more suitable for the numerical implementation within the Lagrangian description framework. The contact law in frictionless contact then reads 8t 2 ½O; T; 8X 2 S ð1Þ 0c ; gðX; tÞ 6 0 : if gðX; tÞ < 0; then T N ðX; tÞ 0 : if gðX; tÞ 0; then T N ðX; tÞ P 0 ð20Þ

In this work, use will be made of the penalty method and the above contact law will be replaced by a smoother one

T N N hgi ð 21Þ
where N > 0 is the penalty parameter and the Macaulay brackets hgi are defined as

hgi g if g P 0 0 if g < 0 ð22Þ
One can check that the contact laws ( 20) and ( 21) are objective by combining Relations ( 13) and ( 18).

The constrained minimization problem for the deployment

Potentials

The internal potential energy of elastic strains P int is defined from the strain energy density wðEÞ via

P int ðUÞ Z S 0 HwðEÞdS 0 ð23Þ
The external loadings are of two types: (i) a follower force which is the internal pressure pðtÞ inflating the membrane struc ture during its deployment and (ii) a dead load such as the weight of the membrane itself. Since the pressure is applied on a closed surface, it is a conservative load like the dead load. Thus, all the external loadings give rise to the potential

P ext ðU; tÞ Z S 0T T Á UdS 0 þ pðtÞVðUÞ ð 24Þ
where T is the dead load per unit surface prescribed over S 0T ; VðUÞ is the volume enclosed by the membrane surface, which depends on the current position of the membrane and can be expressed as a function of the displacement field U by VðUÞ 1 3

Z S x Á ndS 1 3 Z S ðX þ UÞ Á ndS ð25Þ
where n is the outward unit normal vector to the membrane sur face, as shown in Fig. 1. In fact, the membrane considered in this work is very thin so that its weight is almost negligible. One then has to add some other dead loads which will be specified later in the numerical and experimental parts in order to stabilize the quasi static deploy ment process.

From the contact law (21), the potential of the contact forces P c is a quadratic functional of hgi:

P c ðUÞ Z S ð1Þ 0c 1 2 N hgi 2 dS 0 ð26Þ

The constrained minimization problem

To determine the solution of the problem, one may use the well known principle of minimum total potential energy which states that the actual displacement field U minimizes the total potential energy PðU; tÞ P int ðUÞ P ext ðU; tÞ þ P c ðUÞ among all kinematically admissible displacements. The energy minimization method was successfully applied by [START_REF] Bonet | Finite element analysis of air supported membrane structures[END_REF] to air supported membrane structures and [START_REF] Gil | Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains[END_REF] to prestressed membranes undergoing large strains.

Here, the minimization problem will be transformed into an other one which proves to be more convenient for the deployment computations. Indeed, as we shall see in the experimental works (Section 7), the inflating pressure pðtÞ is not a monotonous function during the deployment: it first increases up to a maximum value and then decreases as soon as the deployment is in a sufficiently advanced stage. There is not a one to one relationship between the internal pressure and the maximum displacement for instance. Meanwhile, the volume VðUÞ enclosed by the membrane surface keeps on increasing till the end of the deployment.

As a consequence, one had better take volume VðUÞ as the con trol parameter, instead of pressure pðtÞ. To do this, let us rewrite the total potential energy as

PðUÞ Z S 0 HwðEÞdS 0 Z S 0T T Á UdS 0 þ Z S ð1Þ 0c 1 2 N hgi 2 dS 0 ð27Þ
where only the potential due to external dead loads is retained in P ext and there is no more potential due to the pressure. Solving the problem then reduces to the following minimization problem of PðUÞ under constraint: given a prescribed volume V, find U minimizing PðUÞ

Z S 0 H wðEÞdS 0 Z S 0T T Á UdS 0 þ Z S ð1Þ 0c
1 2

N hgi 2 dS 0 subject to hðUÞ

1 3 Z S x Á ndS V 0 ð28Þ
Besides the volume constraint hðUÞ VðUÞ V 0, there exist other constraints to be taken into account as well, such as the boundary conditions and the symmetry conditions (only half the tube is modeled in finite element computations). In Relation (28), only the volume constraint is mentioned since it is less standard than others and deserves particular attention in the deployment problem. The other constraints are not mentioned explicitly for the sake of clarity.

Algorithm of the projected gradient method

Minimization problems are usually solved by means of either first order methods [START_REF] Gil | Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains[END_REF] or second order methods [START_REF] Bonet | Finite element analysis of air supported membrane structures[END_REF][START_REF] Gil | Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains[END_REF]. Second order methods necessitate the tan gent stiffness matrix which is often singular when the membrane is not yet taut. Although first order methods have only a linear rate of convergence, they prove to be particularly efficient for inflatable structures as they require computing the first gradient of the total potential energy only, not the Hessian matrix. Thus, it is possible to compute equilibrium configurations which exist when the membranes are not taut. Here, use will be made of the projected gradient method formulated by [START_REF] Goldstein | Convex programming in Hilbert space[END_REF] and [START_REF] Levitin | Constrained minimization methods[END_REF] in order to solve the minimization problem under volume control (28). This is a first order method of the same type as those used in [START_REF] Bouzidi | Numerical solution of hyperelastic membranes by energy minimization[END_REF] and [START_REF] Gil | Structural analysis of prestressed Saint Venant-Kirchhoff hyperelastic membranes subjected to moderate strains[END_REF].

Let PðUÞ be the functional to be minimized under the con straints hðUÞ ðh 1 ðUÞ; h 2 ðUÞ; . . . ; h n ðUÞÞ 0. The simple update relation U iþ1 U i þ a s i , where s i is the descent direction at itera tion i and a the descent step, needs to be corrected since it does not fulfill the constraints hðUÞ 0. The key idea of the algorithm lies in using the projected point U i of U i onto the constraints sub space hðUÞ 0 and searching the solution increment U iþ1 U i in the tangent hyperplane of hðUÞ 0, see Fig. 4. Thus, the update relation writes

U iþ1 U i þ as i þ ðU i U i Þ ð 29Þ
where the descent direction s i is chosen to be a unit vector tangent to the constraints subspace.

Therefore, problem (28) with unknown U is transformed into another one whose unknown is the descent direction s:

Find s minimizing s T $ U P subject to B T s 0 and s T s 1 ð30Þ
where B $ U h is the gradient matrix of h. Here, the functional to minimize represents the first variation of energy P defined in (27). The Lagrangian L associated to problem (30) is

Lðs; k; lÞ s T $ U P s T Bk lðs T s 1Þ ð 31Þ

whence the stationarity condition of L @L @s

$ U P Bk 2ls 0 ð32Þ
Premultiplying Relation (32) by B T gives the Lagrange multiplier k:

k B T B h i 1 B T $ U P ð33Þ
Re inserting k in (32) then leads to s:

s 1 2l I B B T B h i 1 B T $ U P ð34Þ
where the multiplier l has to be adjusted so that ksk 1. Further more, it can be shown from the relation

hðU i Þ B T U i U i that U i U i BB T h i 1 BhðU i Þ ð 35Þ
Finally, the update scheme (29) becomes

U iþ1 U i þ a 2l I B B T B h i 1 B T $ U P BB T h i 1 BhðU i Þ ð 36Þ
The unknown a is obtained by minimizing the energy PðU iþ1 Þ PðU i þ as i Þ with respect to scalar variable a by the so called line search method. The iterative process is stopped when the displacement increment kU iþ1 U i k is small enough.

Finite element implementation

The projected gradient method is implemented by means of the finite element method. We shall focus on the discretization of the contact potential (26) and the volume (25), as other terms in (28) can be dealt with in a standard way. In the sequel, the following notational convention is adopted for the matrix representations: curly brackets { } designate a column vector and the square brack ets [ ] a general matrix.

Let the contactor surface S ð1Þ 0c be subdivided into a number of isoparametric finite elements e ð1Þ (the superscript '(1)' is used throughout to remind of the contactor surface S ð1Þ 0c ), then the con tact potential (26) can be recast as the sum of the element terms related to elements e ð1Þ and evaluated using the Gauss quadrature rule:

P c X e ð1Þ N 2 Z e n hgi 2 jðnÞdn 1 dn 2 X e ð1Þ N 2 X Gauss points hgi 2 jðnÞwðnÞ ð37Þ
where e n is the reference element of e ð1Þ ; n ðn 1 ; n 2 Þ 2 R 2 the parent coordinates of point X 2 e ð1Þ (see Fig. 3) and jðnÞ kX; n 1 ðnÞ Â X; n 2 ðnÞk is the area Jacobian. In the last sum in Relation (37) all terms are computed at Gauss point n and wðnÞ is the associated weight.

In any element e ð1Þ & S ð1Þ 0c , the reference and current coordinates are interpolated by fXg ½N ð1Þ ðnÞfXg e ð1Þ fxg ½N ð1Þ ðnÞfxg e ð1Þ ð38Þ

In (38) b for instance, fxg is the column vector containing the three coordinates of point x in a fixed Cartesian orthonormal basis ðe 1 ; e 2 ; e 3 Þ. The column vector fxg e ð1Þ contains the nodal coordinates of element e ð1Þ , and the matrix ½N ð1Þ ðnÞ contains the shape func tions of element e ð1Þ .

In a similar way, the target surface S

0c is subdivided into a num ber of isoparametric finite elements. Given a quadrature point X in element e ð1Þ , i.e. x /ðX; tÞ, let e ð2Þ the (or an) element of S x. Both element e ð2Þ and its node number, which depend on the projection y, vary as functions of the quadrature point X. Moreover, one element e ð1Þ can be in contact with several elements e ð2Þ ; consequently, the whole set of quadrature points of a given element e ð1Þ in (37) generally involves several elements e ð2Þ in S ð2Þ 0c . It results that the approach presented herein is of the node to segment rather than the segment to segment type (for the mentioned types, see, e.g., [START_REF] Kikuchi | Penalty finite element approximation of a class of unilateral problems in linear elasticity[END_REF][START_REF] Harnau | Algorithmic aspects in large deformation contact analysis using 'solid-shell' elements[END_REF][START_REF] Puso | A mortar segment-to-segment frictional contact method for large deformations[END_REF][START_REF] Puso | A mortar segment-to-segment frictional contact method for large deformations[END_REF].

By using Relations (38) and ( 39), the potential of the contact forces ( 26) is computed from where e n is the reference element of e, the current position x is interpolated as in (38) by fxg ½NðnÞfxg e and a a ½N ;a ðnÞfxg e .

P c X e ð1Þ

Contact detection

The main task in evaluating the contact potential (40) and its gradient is to determine the projection y of an integration point x onto the target surface S ð2Þ c . The contact detection is a nonlinear procedure which may be highly CPU time consuming and requires an iterative process as explained in [START_REF] Benson | A single surface contact algorithm for the postbuckling analysis of shell structures[END_REF]. In order to reduce the computational cost, most algorithms involve a preliminary vicinity test, which consists in identifying the closest elements likely to interpenetrate so that the contact computation can be restricted to those elements only, for more details see [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF].

It should be mentioned here another difficulty specific to the deployment of membrane structures. Usually, the contactor sur face S ð1Þ c may come into contact with the target surface S ð2Þ c from one side of the target body, namely the outward side as defined by the outward normal vector n ð2Þ to S ð2Þ c , see the convention about normal n ð2Þ in definition ( 9). This situation is shown in Fig. 5a for the points x; y and the normal m in Fig. 3. On the other hand, when the contact involves hollow bodies as is the case with (closed) membranes, the target surface S ð2Þ c may encounter the contactor surface S ð1Þ c from either side [START_REF] Wang | The inside-outside contact search algorithm for finite element analysis[END_REF]. As can be seen in Fig. 5b, point y may potentially come into contact with point x or x 0 , located in both sides of S ð2Þ c with respect to normal m.

Point x is in the same situation as in Fig. 5a and is treated as usual. Point x 0 needs a special operation: in order to prevent it from cross ing point y, one has to replace normal m by the opposite normal m.

Validation of the contact algorithm

The contact algorithm is tested on the following simple exam ple which admits of an analytical expression for the contact sur face. Consider two identical spherical balloons of initial radius R inflated by an internal pressure p and compressed by two opposite forces of same magnitude F passing through the centers of the spheres, as shown in Fig. 6a. The problem is symmetrical with re spect to the mid plane perpendicular to forces F and the contact surface is a planar disk. The balloons are modeled as membranes (thus having no bending stiffness) and the contact is assumed to be frictionless. By writing the equilibrium balance for one of the two balloons with its air and without the surface of contact, see Fig. 6b, one finds that the contact area S c is simply the ratio of the compressive load to the internal pressure: S c F=p, thus inde pendent of the geometry and the constitutive material.

The numerical computations are done with the radius R 50 cm, the membrane thickness H 50 lm and by assuming the Saint Venant Kirchhoff potential wðEÞ lE

: E þ k 2 tr E ð Þ 2 , so that it comes from (6): R 2lE þ ktrðEÞI () E 1 þ m E R m E trðRÞI ð42Þ
where I is the unit second order tensor, k and l are the Lamé con stants, E is the Young's modulus and m the Poisson's ratio. Here, we take E 5:10 7 Pa; m 0:3. For symmetry reasons, only a quarter of the system is discretized, the mesh is made of three node triangle elements as shown in Fig. 7. The penalty parameter is and internal pressure p increasing from 350 Pa to 850 Pa. At each computational step, we search the whole surface of the balloons and as soon as an element is found to be in contact, its entire area is added to the contact area. Thus, the numerical value of the con tact area can only vary in a discontinuous way, however fine the mesh is and however small the computational step is. It is difficult to precisely evaluate the contact area and the mesh chosen here is fine enough for the discontinuity to be acceptable. The results are displayed in Tables 1 and2 for the numerical contact area in com parison with the theoretical values, together with the axial dis placement u of the point of application of force F. In the second calculation, the contact area decreases as pressure p increases, it is the volume of the two balloons that increases with p.

Because of the above mentioned discontinuity, the error in creases when few finite elements are in contact, i.e. when the load is small or the pressure is high. The numerical values are in good agreement with the analytical ones: the maximum error is less than 4% in the first calculation and 3% in the second one. This example shows that the numerical developments for deployment by contact and self contact work well.

Numerical computations of the deployment

The numerical computations of the quasi static deployment are performed on a folded cylindrical membrane tube (length = 80 cm, diameter = 15 cm or 20 cm, thickness = 50 lm), made of a hyperelastic St Venant Kirchhoff material with Young's modulus E 5:10 9 Pa and Poisson's ratio m 0:3.

Initially, the tube is completely flattened and then folded in two equal parts as shown in Fig. 8. The bottom edge is fixed while the top edge is constrained to move freely along the z axis perpendic ular to the flattened position of the tube. At each step of the deployment, the volume of the tube is set to a prescribed value as stated in (28), which amounts to inflate the tube with a certain internal pressure. The computations are carried out in quasi statics and a dead load Fz is uniformly distributed on the top edge in or der to stabilize the deployment process. Here, the dead load is due to a mass of 691.7 g, 1383.9 g, 2050.7 g or 2742.9 g, i.e. F 6:78 N;13:62 N;20:11 N or 26.9 N.

Because of the symmetry of the problem with respect to the middle vertical plane, only half of the tube is modeled. A very small initial gap of 0.01 cm is prescribed between the top and the bottom edges (equivalent to an opening angle of about 0.01°), so that the contact surfaces are very close yet not in contact at the first com putational step. Fig. 8 shows the mesh used containing 529 nodes and 1024 three node triangular elements. The penalty parameter for the contact is N 10 6 N=m 3 in all the numerical computations.

As noticed by [START_REF] Salama | Simulation of the deployment dynamics of inflatable structures[END_REF] in their study of the dynamic deployment of inflatable structures, high values of the Young's modulus may cause a numerical instability. This kind of difficulty is commonly observed in the deployment computations. In the present work, one finds a similar numerical locking phenomenon, which can be accounted for by two main reasons: (i) the facettiza tion of the membrane surface due to the finite element mesh, and (ii) the significant difference between the in plane and out of plane stiffnesses inherent to the membrane model. With a high value of the Young's modulus E, the strain energy is found to be prevalent compared to other energies and generates a numerical noise, leading to a premature convergence towards local minima of the total potential energy. This noise may decrease with the mesh refinement but it remains a major hindrance to the convergence. To remedy the situation, a simple solution is to arti ficially and temporarily diminish the membrane stiffness of the tube, by reducing E in the first iterations of each computational step and resetting it to its real value in the final iterations of the step. The choice of the artificial value for E has little influence on the deployment kinematics since the deployment involves large displacements and rotations, rather than large strains. For our example considered herein, we set E 5:10 7 Pa in the first itera tions a small value conveniently chosen so as to avoid the above mentioned numerical difficulties and yet maintain small strains and reset the Young's modulus to its real value E 5:10 9 Pa in the final iterations of the step.

Fig. 9 shows the 20 cm diameter tube at different stages of deployment with applied dead load F 13:62 N. The numerical re sults show that the deployment is possible only if the diameter of the tube is large enough or load F is small enough: for the heaviest load F 26:9 N (2742.9 g), convergence is achieved for the 20 cm diameter tube, not for the smaller tube. The inflating pressure p is plotted versus the top displacement along direction z in Fig. 10 for the prescribed F values. As can be seen, for a given dead load F, the pressure reaches a maximum value at about halfway and then decreases until the end of the deployment. The higher the load F, the higher the maximum pressure is. On the other hand, the vol ume of the tube is monotonically increasing throughout the deployment process.

Deployment of a Z folded tube

Another example is briefly described here to show how the numerical computations work on more complex geometries with multiple self contact surfaces. We consider the same tube as in the previous section and assume now that it is folded in four equal parts. As before, only half of the tube with respect to the middle vertical plane is modeled. The boundary conditions are prescribed so that the top and medium edges of the tube slide freely along the vertical direction and a dead load force F 1 N is applied on the top edge. The finite element model is shown in Fig. 11.

As may be seen in Fig. 12, the upper and lower parts of the tube are deployed in a symmetrical way in the first stages of the deploy ment. However, the symmetry is broken as soon as one crosses the limit point corresponding to the maximum pressure on the deploy ment curve. The displacement field of the tube then bifurcates from the symmetrical branch and the lower half of the tube is com pletely deployed before the upper one. Further investigations could be made to control the symmetrical and asymmetrical branches in the bifurcation context but they fall out of the scope of this work.

Experimental validation

This section describes the experimental works conducted in order to validate the numerical results obtained in the previous Section. Experiments on deployment have been investigated by a large number of authors in the dynamical context, see, e.g., Free land and Veal (1998), and [START_REF] Miyasaki | Deployment dynamics of inflatable tube[END_REF]. The study presented here differs from these works in that the deployment is conducted in a quasi static way, as done in [START_REF] Buytet | Deployment of membranous tubes by air inflation at low pressure[END_REF]. We consider the membrane tube folded in two, having the same geometries and material properties as in the numerical part. The tube has a length of 80 cm, a diameter of 15 cm or 20 cm, it is made of a flexible and isotropic thin film (50 lm thick) with Young's modulus E 5:10 9 Pa and Poisson's ratio m 0:3. At the beginning, the tube is completely flattened and folded in two. During the deployment, the bottom of the tube is kept in fixed position by a pivot pin and the top is attached to an airborne car riage mounted on a fixed vertical axis, see Figs. 13 and 14. The air borne carriage is designed so as to allow the top of the tube to move freely with as little friction as possible. A vertical dead weight F is applied on the top of the tube in order to stabilize the deployment in quasi statics.

The deployment takes place as air is injected in the tube and the tube increases in volume. For each given quantity of air, the posi tion of the top is measured when the internal pressure inside the tube is stable. The quantity of the inflating air is increased until the tube is almost completely deployed, namely when the top dis placement is about 65 cm. From this moment on, the tube is grad ually deflated until it returns to the initial flattened and folded configuration. Fig. 15 displays the pressure versus the top displace ment for the 15 cm and 20 cm diameter tubes, under four dead loads F = 6.78 N, 13.62 N, 20.11 N and 26.9 N. For a given geometry and dead load, several tests are carried out, giving rise to slightly different response curves in Fig. 15.

It is noteworthy that wrinkles may appear virtually randomly all over the membranes throughout the deployment process. For a given geometry and dead load, the distribution of the wrinkles changes from one test to the next; however, this has insignificant influence on the global response of the deployment.

As can be seen in Fig. 15, the inflation and deflation curves do not coincide, giving rise to a hysteresis, contrary to what can be ex pected for a conservative system. For a given position of the tube top, the pressure in the inflation stage is always higher than that in the deflation stage. The same phenomenon is observed in all tubes, for all applied dead loads and for various materials used.

One may think that the hysteresis is due to the fact that some friction can exist between the airborne carriage and the vertical axis, or over the surfaces where the membrane comes into self contact. However, that should not be the cause since a lot of pre cautions have been taken in the experimental settings in order to diminish the friction: in addition to the above mentioned air cush ion system, the surfaces of the membrane are lubricated with oil. Also, the vertical dead load applied to the top of the tube is pre cisely measured by means of a force sensor inserted between the dead load and the top of the tube (see the load gauge in Fig. 13).

The measured load is plotted in Fig. 16 in the case of the 20 cm diameter tube submitted to F = 6.78 N and 26.9 N. It is found that the measured load is equal to the theoretical load F in the inflation stage, but it is slightly lowered in the deflation stage, about 0.4 N for F = 26.9 N. Additional masses are then used to compensate for this difference so that the same dead load is actually applied dur ing the inflation and the deflation, as shown in Fig. 17. In spite of this correction, the hysteresis observed in Fig. 15 remains un changed. The conclusion is that the friction introduced by the car riage, if any, cannot be accountable for this hysteresis. A close look at the membrane surface reveals that the wrinkles are pronounced in the beginning of the inflation stage and most of them disappear when pressure p gets across its peak value. In the subsequent deflation stage, the membrane is always smoother with fewer wrinkles. The wrinkles appearing in the first stage are so marked that the membrane could be plastified at these locations and the mechanical energy used to smooth out the wrinkles may be the cause of the observed hysteresis.

Fig. 18 shows that the numerical results obtained in Section 6 agree very well with the experimental results. In particular, the limit point in the inflating pressure p encountered at halfway of the full deployment is predicted with a good precision. The differ ence between the numerical and the experimental pressure in creases with the dead load; its maximal value is about 10%.

Conclusions

The deployment of inflatable structures in quasi static condi tions is a difficult problem as it involves different types of nonlinear ities such as the large deformations, the inflating pressure which is a follower loading, the contact and self contact phenomenon; and most importantly the difficulty due to the singularity of the struc tural stiffness during the deployment stage. The minimization formulation combined with the method of projected gradient has been proposed as a simple tool to solve the problems of this type. The efficiency of the approach has been proved through the case of a simply folded tube as the numerical results obtained have been found to be in good accordance with the experimental results.

  Fig. 1. Parametrization of the membrane middle surface. Fig. 2. Parametrizations of the contact surfaces.

S

  ðiÞ c , whereas the notation x designates a particular point of interest on S ð1Þ c and y the opposite point of x on S ð2Þ c defined by (10). The point x 2 S ð1Þ c being associated to (at least) one point y 2 S ð2Þ c , the proximity g is defined by g m:ðx yÞ ð 11Þ

  h ð1Þ ðnÞ#x /ðX; tÞ j ð1Þ ðn; tÞ ð14Þ whereas points g; Y and y are linked by A region of R 2 ! S ð2Þ 0c ! S ð2Þ c g#Y h ð2Þ ðgÞ#y /ðY; tÞ j ð2Þ ðg; tÞ ð15Þ By making x ð2Þ j ð2Þ ðn ð2Þ ; tÞ in (10), one finds that g is character ized by g argmin n ð2Þ kx j ð2Þ ðn ð2Þ ; tÞk

  reference area, dS ð1Þ its spatial counterpart), the nominal traction T can be resolved in terms of normal m at y 2 S ð2Þ c as follows TðX; tÞ T N ðX; tÞm ð19Þ where the material and spatial normal traction components are related by t N dS ð1Þ T N dS ð1Þ o .

Fig. 3 .

 3 Fig.3. A contactor point x and its the projection y on the target surface. In comparison with Fig.2,x is a particular point x ð1Þ and y a particular point x ð2Þ .

  point Y corresponding to the projection y proj S ð2Þ c x 2 S ð2Þ c . In element e ð2Þ , the reference and current coordi nates are interpolated by fYg ½N ð2Þ ðgÞfYg e ð2Þ fyg ½N ð2Þ ðgÞfyg e ð2Þ ð39Þ where g ðg 1 ; g 2 Þ 2 R 2 is the parent coordinate of Y and the super script '(2)' is used to remind of element e ð2Þ & S ð2Þ 0c . The column vec tor fyg e ð2Þ contains the nodal coordinates of element e ð2Þ and the matrix ½N ð2Þ ðgÞ contains the shape functions of element e ð2Þ . Given a quadrature point X in element e ð1Þ & S ð1Þ 0c , there is at least one element e ð2Þ & S ð2Þ 0c containing the reference point Y of y proj S ð2Þ c

ð40Þ

  The reference middle surface S 0 is discretized into isoparametric elements e. The volume (25) of the structure is 1 Â a 2 Þdn 1 dn 2 ð41Þ4. Projected gradient algorithm.

N 10 6 Fig. 5 .FFig. 7 .

 657 Fig. 5. (a) Usual contact condition; (b) contact from either sides in deployment.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Deployment of a folded membrane tube.

Fig. 10 .

 10 Fig. 10. Internal pressure versus top displacement.

Fig. 12 .

 12 Fig. 12. Deployment stages of the Z-folded tube with applied dead load F 1 N.

Fig. 13 .

 13 Fig. 13. Experimental schema.

Fig. 16 .Fig. 17 .

 1617 Fig. 15. Experimental deployment curves.

Table 1

 1 Analytical and numerical contact areas versus the applied load (pressure p 1000 Pa). Analytical and numerical contact areas versus the internal pressure (force F 10 N).

	Force	Contact area Scðm 2 Þ		Top
	F (N)	Analytical values	Numerical values	Difference (%)	displacement u (m)
	20	0.020	0.0203	1.5	0.1041
	30	0.030	0.0310	3.6	0.1753
	60	0.060	0.0613	2.2	0.3814
	70	0.070	0.0693	1.0	0.4571
	Table 2
	Pressure p	Contact area Scðm 2 Þ		Top displacement
	(Pa)				u (m)
		Analytical	Numerical	Difference
		values	values	(%)
	350	0.0286	0.0290	1.6	0.171
	400	0.0250	0.0256	2.5	0.151
	450	0.0220	0.0224	2.0	0.134
	500	0.0200	0.0202	1.3	0.120
	550	0.0182	0.0182	0.1	0.107
	650	0.0154	0.0158	2.7	0.085
	750	0.0133	0.0137	2.7	0.066

Acknowledgement

The National Space Center (Centre National d'Etudes Spatiales, CNES) in Toulouse, France, is gratefully acknowledged for having partly supported this work. The authors would also like to thank K. A. Brakke for his Surface Evolver freeware (Brakke, 1992) which has been very helpful throughout all numerical developments and implementations in this paper.


Alexandre Piaget

Matthieu Museau

Henri Paris

Manufacturing Space Homogeneity in Additive Manufacturing -Electron Beam Melting Case

Keywords: Quality Management, Additive Manufacturing, Electron Beam Melting, Manufacturing Space Homogeneity

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The additive manufacturing technologies have become essential tools for modern industry [START_REF] Gibson | Development of Additive Manufacturing Technology[END_REF] [START_REF] Gao | The Status, Challenges, and Future of Additive Manufacturing in Engineering[END_REF]. Among the different additive manufacturing technologies, EBM (Electron Beam Melting) is able to manufacture a wide range of metallic parts (massive, topologically optimized, lattice structure) [START_REF] Hrabe | Effects of Processing on Microstructure and Mechanical Properties of a Titanium Alloy (Ti-6Al-4V) Fabricated Using Electron Beam Melting (EBM)[END_REF] [START_REF] Novak | Auxetic Cellular Meterialsa Review[END_REF]. Compared to SLM (Selective Laser Melting), the EBM technology is still poorly studied [START_REF] Sing | Laser and Electron-Beam Powder-Bed Additive Manufacturing of Metallic Implants: A Review on Processes, Materials and Designs[END_REF]. However, there is a strong need to further master this technology. Aeronautic, spatial and also medical fields are particularly interested in mastering the quality provided by this process [START_REF] Daneshmand | Aghanajafi: Description and Modeling of the Additive Manufacturing Technology for Aerodynamic Coefficients Measurement[END_REF] [START_REF] Waller | Nondestructive Evaluation of Additive Manufacturing[END_REF][8] [START_REF] Petrovic | Additive Manufacturing Solutions for Improved Medical Implants[END_REF].

In order to control the quality, without destruction of manufactured parts, test specimens are manufactured simultaneously [START_REF] Mat -Str | PRINTtemps de la Fabrication Additive[END_REF]. By testing those specimens, information about the quality of manufactured parts is obtained. Most of the time, specimens are arranged in the edges of the manufacturing space. A question arises about the homogeneity of the manufacturing space. Indeed, heterogeneity could lead to errors in the manufacturing and control processes. This key point that is not addressed in the literature, is the topic of this paper.

The problem is approached by studying the geometrical deformation which is the main defect. Previously, several parts showed defects when they were manufactured within certain areas of the manufacturing space. Thus, there is a need to characterise this space and its heterogeneity, and the way it impacts on the manufactured parts. These characteristics are related to the machine used (described below), but the experimental method can be applied to other machines.

In the next part, the context of the study and the experimentation are exposed. The type of deformation and measures are presented. Secondly, the possible sources of this defect are analyzed. We also show the impact they may have on manufactured parts. Finally we conclude about the dangers of this phenomenon and the prospects for resolution.

Defect identification

Context

The Arcam AB A1 machine [START_REF] Arcam | Creating New Opportunitites in Design and Production[END_REF] [START_REF] Arcam | Arcam A1, the Future in Implant Manufacturing[END_REF] uses an electron beam to melt powder (Figure 1). The powder from the hoppers is spread over a platform in 50µm thick layers with a rake. The electron beam sinters the whole of the introduced powder. This increases the thermal and electrical conductivity. Then it locally melts the powder to manufacture the desired geometry. The platform finally goes down of the thickness of a layer and this cycle can be repeated.

Figure 1: Arcam EBM systems, schematic architecture [START_REF] Arcam | Creating New Opportunitites in Design and Production[END_REF] In order to qualify the geometric defect, some parameters have to be chosen. Ti-6Al-4V powder from Arcam AB is used for the study. Vayre [START_REF] Vayre | Conception pour la Fabrication Additive, Application à la Technologie EBM[END_REF] shown that Arcam AB standard parameter set is well balanced for printing several kind of geometry (lattice structures, massive workpieces…) at the same time. Thereby, this parameter set1 is kept for the whole experimentation. Otherwise, workpieces will be built without any supporting strategy. Indeed, EBM technology doesn't necessarily involve a supporting strategy to build workpieces [START_REF] Obaton | Fabrication additive: Etat de l'Art et Besoins Métrologiques Engendrés[END_REF]. Based on Vayre [START_REF] Vayre | Identification on Some Design Key Parameters for Additive Manufacturing: Application on Electron Beam Melting[END_REF], there is no need to add support to the experimental workpieces (described in the section below) because of the experimented workpieces thinness. Hereafter, evidence that supports aren't necessary to manufacture these workpieces is provided.

For the experimentation, the largest manufacturing space available with the machine is needed. But the manufacturing space is limited by the use of software2 that reduce the initial manufacturing space (200*200*180) to 195mm in length, 195 in width and 180 in height.

Experimentation

The first objective of this study is to evaluate the impact of the position in the manufacturing space on the final workpiece geometry. Indeed, deformations are observed on several workpieces made by EBM. In order to show and identify this phenomenon, the following experimentation is realized.

The principle is to build one workpiece at different locations in the manufacturing space. Then, differences and defects can be observed and correlated to the location. Two kinds of workpieces have been selected for the test (Figure 2). The first one is a massive workpiece. The second one is a lattice structure. Thus, those workpieces represent the two kinds of geometry cases encountered. The height and the diameter of the pipes are 30mm. The thickness of the wall for the massive pipe and of the beam for the lattice one is 1mm. Inside each pipe, two other pipes with different diameters (Ø20mm and Ø10mm) are integrated to the experimentation in order to have different thermal conditions. Indeed, the Ø10mm workpiece always has a Ø20mm and a Ø30mm pipes for neighbours while the Ø30mm workpiece has a number of neighbours depending on its position.

The experimentation consists in testing 25 locations of the manufacturing space. Grouped by three, workpieces will be repeated on each locations of a unique altitude of the manufacturing space in order to avoid scan length problems [START_REF] Safdar | A study on Electron Beam Melted Ti-6Al-4V[END_REF]. By measuring each workpieces, defects are characterized and linked to the location.

Observations show that the defect appears in first layers. Therefore, the lower part of the workpieces is analysed and that is the reason why the workpieces are represented upside down later in the document. The measurement protocol for the massive workpiece uses a three-dimensional optical control machine 3 in order to measure deformation on the first manufactured layer. This layer is controlled at 24 points (every 15 degrees) distributed all along the surface. The result of a measure is the height of the controlled point with a 2µm precision. From this, we recovered the difference between the nominal height and the measured one (Figure 3, left). Measuring the first layer deformation does not provide information on the deformation of the following layers. With the lattice structure, the target is to measure the defect spread through several layers. Consequently, the measurement protocol has to be adapted. A tool has been created to control the external geometry of the lattice pipe (Figure 3, right). This way, we may have information about the defect on the pipe: its angular position and altitude.

Results

Measures are presented in a radar chart in order to represent the size and the angular position of the defect. So, in the diagrams, the more a point is far from the middle, the more the defect is important. As shown in Figure 4, defect occurs on some workpieces (on both pictures, the left workpiece fits the CAD model and the right includes a deviation from the model). The defect is materialized in two forms: material loss and geometric deformation. Lattice structures show both failure (Figure 4): between the well-shaped and the defective workpiece, deformations are observed in first layers (missing) and in diamond meshes (misshapen).

Figure 4: massive (left) and lattice (right) workpieces placed upside down Figure 5 presents the results for the massive workpiece. From these results, occurrences of defects on nearby workpieces of borders are observed. Furthermore, the defect is more critical as the point observed on a workpiece is close to the boarders. The biggest gap measured between the nominal and the manufactured surface is 3mm high. Whereas, for the inner workpieces, measured defects are smaller than 0.8mm. But this fact corresponds with the 50µm roughness owned by the material resource [START_REF] Safdar | A study on Electron Beam Melted Ti-6Al-4V[END_REF]. Figure 6 presents results for the lattice workpieces in numbers of defective cells, for an angular position and starting from the bottom of a workpiece. The same conclusion can be drawn as there is no defect in the inner workpieces and the external workpieces have an outwardly directed defect. In addition, the highest measured defect is 8mm high from the first layer. Another interesting point is that a group of workpieces (located at the bottom right in Figure 6) has undergone a 2mm mechanical shifting. The 3 curves are overlapping because workpieces cells are all affected in the same way. What is interesting in this point is that it brings other elements to explain and understand the sources and consequences of the deformation phenomenon, although this shifting is not repeatable. Figure 7 is a photography that shows the encountered problem. This shifting is in the movement direction of the rake. Finally, observations and measures show that only the workpieces near the boarders of the manufacturing space are reached by the defects. In corners, the defect reaches its highest value: the maximal height of the defect is 8mm from the first layer. Those defects are significant enough to justify a thorough study with the aim to master and avoid this problem.

Phenomenon Analysis

Hypothesis

To explain the phenomenon, several hypotheses have been expressed. First, the electron beam hypothesis has been selected. Indeed, to reach boarders, the electron beam has to undergo an important deflection that leads to change the shape of the beam from a circle to an ellipse. The energy broadcasted to the powder is lower as the surface of transmission gets bigger. It generates a decrease in the temperature that prevents the melting of the powder. This decrease might explain the location of the defect and the defect type also.

The second hypothesis concerns the workspace in the machine. At the borders of the manufacturing space, next to the sintered powder, there are cold powder and the metal enclosure that might absorb the energy of the sintered powder and thereby generate a loss of temperature. As TA6V is a bad thermal conductor [START_REF] Budzik | Analysis of the Accuracy of Reconstructed Two Teeth Models Manufactured Using the 3DP and FDM Technologies[END_REF], this phenomenon might just take place at the border between the sintered powder and the cold powder. Lowering temperature locally leads to the same conclusion as the electron beam hypothesis.

A third hypothesis arises from the mechanical shifting (Figure 7). As a portion of the workpiece has been moved in the rake direction, this hypothesis involves the rake: when the rake spreads the powder out, it can move the matter underneath by the thrust applied on the powders. That might cause deformation and matter loss in the workpiece.

The first and the second hypothesis might explain the difference between the workpieces but they don't solve the matter shifting. On its side, the third hypothesis elucidates the shifting but there is no reason for the nonhomogeneity of that phenomenon.

The hypothesis we kept is a combination of the previous ones: the contact area of the electron beam is circular in the powder bed center, but away from the center, it becomes elliptical and provides less energy to the powder. In addition, the energy leaks with the proximity to cold powders and metal enclosure. This makes the sintering of the powder weaker near the boarders. Then, the rake passage might easily move the poorly sintered powder and the molten material suspended in the powder. This might lead to a loss of material as well as deformation for the manufactured parts. In addition, the sintering weakness modifies the thermal conductivity, the cooling of the manufactured parts is no more homogeneous and it generates additional deformation.

Analysis

In order to validate the hypothesis, an analysis of the sintered powder is made. This analyse was conducted with an X-ray microtomography machine. The machine is able to reconstruct a 3D model of the studied sample with voxels. This reconstruction aims to display and measure the size of the necks binding in the sintered powder. A difference of necks size between the powders in the centre and in the periphery of the manufacturing space would indicate a significant difference of the sintering [START_REF] Olmos | Etude du Frittage de Poudres par Microtomographie in situ et Modélisation Discrète[END_REF].

The samples of powders are one centimetre side cube. This size allows observation of a large number of powder grain that limits the impact of defects in the powder. Due to the titanium opacity to X-ray, getting a 1µm resolution with titanium parts requires more power than supplied by the machine [START_REF] Thiery | Tomographie à Rayons X[END_REF]. Resolution has been reduced to 5µm in order to provide a clear tomography of the samples. Despite a clear tomography, it becomes impossible to measure necks because of the resolution. However, [START_REF] Olmos | Etude du Frittage de Poudres par Microtomographie in situ et Modélisation Discrète[END_REF] show there is a link between density and sintering progress with spherical powders. Thereby, measuring density permits comparison of the two samples sintering. 1 presents the results of the analysis. With more than 3% difference, assumptions about the difference of sintering are justified. Indeed, according to OLMOS [START_REF] Olmos | Etude du Frittage de Poudres par Microtomographie in situ et Modélisation Discrète[END_REF], 3% are significant enough to make a difference of sintering. This difference impacts the powder stiffness and its heat conductivity. This generates a harmful heterogeneity in the manufacturing space that might have bad effects on manufactured parts [START_REF] Zou | Effect of Sintering Neck on Compressive Mechanical Properties of Porous Titanium[END_REF].

Conclusion

The observed defect has an impact on geometry of manufactured parts with EBM technology. According to their position, parts may undergo geometric deformation and matter loss. It can be inferred from this study that the manufacturing space is heterogeneous.

In the worst case, the 8 first millimeters are affected by this defect. The closest defect to the center is 86mm distant from the center. The more the defect is far from the center, the more important it is. With this information, a first strategy can be used to bypass the defect. Indeed, creating supports 8mm (or more) long all along the first layer permits to transfer the defect from the part to supports. A manufacturing space reduction to an 86mm radius cylinder can also be a solution to avoid the defect apparition. Those 2 strategies have been validated with the workpieces.

As well, modification could be made to change the machine operation and more precisely the electron beam. By bringing more energy, the powder sintering would be homogenized. Improve the insulation of the powder bath with the machine will also limit the energy loss.

Figure 2 :

 2 Figure 2: Experimentation workpieces, massive (left) and lattice structure (right).

Figure 3 :

 3 Figure 3: Representation of the measurement processes: massive (left) and lattice (right) workpieces

Figure 5 :

 5 Figure 5: Example and defect cartography of the massive workpieces (blue: Ø30mm, orange: Ø20mm).

Figure 6 :

 6 Figure 6: Example and defect cartography of the lattice workpieces (blue: 30mm, orange: 20mm, grey: 10mm).

Figure 7 :

 7 Figure 7: Photography of the 2mm mechanical shifting (Ø30mm on left, Ø20mm on right).

Table 1 :

 1 Porosity measured from the tomography reconstruction

	Samples	Centre Periphery
	Porosity (%)	37.84	41.29
	Standard Deviation (%)	1.36	1.92
	Table

Arcam theme melt 50µm for TA6V

Magics (Materialise), Build Assembler (Arcam AB)

Used machine : Vertex, https://www.microvu.com/


Wouter J T Bos
email: wouter.bos@ec-lyon.fr

Robert Rubinstein

Mixing in modulated turbulence. Analytical results

Keywords: Mixing rate, modulated turbulence, passive scalar, isotropic turbulence

Recent numerical results show that if a scalar is mixed by periodically forced turbulence, the average mixing rate is directly affected for forcing frequencies small compared to the integral turbulence frequency. We elucidate this by an analytical study using simple turbulence models for spectral transfer. Adding a large amplitude modulation to the forcing of the velocity field enhances the energy transfer and simultaneously diminishes the scalar transfer. Adding a modulation to a random stirring protocol will thus negatively influence the mixing rate. We further derive the asymptotic behaviour of the response of the passive scalar quantities in the same flow for low and high forcing frequencies.

Introduction

Recent results from direct numerical simulation suggest that the timeaveraged mixing in turbulent flow can be affected by adding a periodic modulation to the energy injection of the turbulent velocity field [START_REF] Yang | Mixing in modulated turbulence: numerical results[END_REF]. Indeed it is shown that the scalar flux towards small scales, or mixing rate, is attenuated by the modulation. On the contrary an enhanced transfer of kinetic energy towards small scales is observed. In the present work we will explain these observations analytically using a second order analysis of simple turbulence closure models.

Closures and turbulence

Nowadays numerical simulations allow us to reproduce academic experiments of isotropic turbulence at Reynolds numbers of the same order as obtained in large wind-tunnel facilities [START_REF] Ishihara | Study of high Reynolds number isotropic turbulence by Direct Numerical Simulation[END_REF]. Closures have therefore partially lost their interest as tools to predict the behavior of turbulent flows at large Reynolds numbers. However, they are still very valuable for the interpretation of results on turbulent flows. If we know the minimum ingredients of a model which is capable of reproducing particular features of a turbulent flow, there is a good chance that we can obtain a qualitative understanding of the particular flow features. If we adopt this viewpoint, it is important to choose the adequate, i.e., the simplest possible model which is able to reproduce the physics of a problem.

The simplest model for spectral transfer of turbulent kinetic energy is probably Kovaznay's closure [START_REF] Kovaznay | Spectrum of locally isotropic turbulence[END_REF]. It should be used if the only features one needs to retain of the nonlinear transfer is the constant flux dynamics of a long inertial range, consistent with Kolmogorov's ideas. Another simple, but more elegant closure was proposed by Leith [START_REF] Leith | Diffusion approximation to inertial energy transfer in isotropic turbulence[END_REF]. In addition to giving the Kolmogorov scaling, his diffusion-approximation is compatible with the thermal equilibrium properties of Galerkin truncated inviscid turbulent flow.

In reference [START_REF] Rubinstein | Reassessment of the classical closures for scalar turbulence[END_REF] the simplest turbulence models were assessed for the case of scalar mixing in isotropic turbulence. Indeed, also in the scalar case, Kovaznay's and Leith's models [START_REF] Leith | Diffusion approximation for turbulent scalar fields[END_REF] are examples of simple models consistent with a conserved flux and local interactions. In the present work we will use such simple closure models to explain the observations of the recent numerical study [START_REF] Yang | Mixing in modulated turbulence: numerical results[END_REF] on mixing in periodically forced turbulence.

More precisely, the questions that we will answer in this manuscript are the following:

• What is the frequency response of the scalar field at low and high frequencies of the modulation of the velocity forcing?

• How is the mixing rate affected when the amplitude of the periodic part of the forcing is not small? Here and in the following, a low (high) frequency refers to a frequency ω of the periodic forcing, small (large) compared to the integral frequency of the turbulent flow.

Notation

In the present work we consider both one-point and two-point statistics, as well as time-, phase-and ensemble-averages, and therefore a word on the notation is helpful in order to clearly distinguish the different types of statistics.

To start with, all quantities we consider are averaged, and since we consider only statistically homogeneous flows, all quantities are independent of the space coordinate. If no tildes or overlines are used, the quantities are ensemble averages, such as k = k(ω, t), depending in principle on the frequency ω and time t. We decompose the ensemble averages in a time-averaged part, indicated by overlined symbols, such as k = k(ω) and the remaining, periodic quantities indicated by a tilde, k = k(ω) (see for instance reference [START_REF] Hussain | The mechanics of an organized wave in turbulent shear flow[END_REF]).

Lower-case symbols (k, ǫ, p) are one-point statistics, upper-case symbols are wavenumber spectra (E = E(κ, ω, t), E = E(ω, κ), ..., with κ the wavenumber). Phase-shifts are indicated by φ = φ(ω), for the one-point statistics and Φ = Φ(κ, ω) for the wavenumber spectra.

Periodically forced turbulence and mixing

One-point statistics

We investigate how turbulent mixing is affected if the velocity statistics contain a periodic component. The energy balance is

∂ t k = p -ǫ, (1)
where k, p, ǫ are ensemble-averages of the kinetic energy, production and dissipation, respectively. The periodic component of the velocity field is created by a modulated forcing,

p = p + p cos(ωt). (2)
The time-averaged and periodic components satisfy,

0 = p -ǫ, (3)
and

ω k sin(ωt + φ k ) = ǫ cos(ωt + φ ǫ ) -p cos(ωt). (4)
In this expression φ k and φ ǫ are phase-shifts. We have assumed that the periodic response of the velocity statistics will be dominated by a contribution at the forcing frequency ω. The frequency responses of k and ǫ were investigated theoretically in [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF][START_REF] Der Heydt | Response maxima in modulated turbulence[END_REF], numerically in [START_REF] Kuczaj | Response maxima in timemodulated turbulence: Direct numerical simulations[END_REF][START_REF] Kuczaj | Turbulence modification by periodically modulated scale-dependent forcing[END_REF] and experimentally in [START_REF] Cadot | Observation of resonances in modulated turbulence[END_REF].

The scalar is forced by an injection p θ = p θ , without periodic component.

But since the velocity field which advects the scalar field is modulated, a trace of the modulation can be found back in the scalar statistics. The scalar variance evolves as

∂ t k θ = p θ -ǫ θ . (5)
The mean and periodic parts of the scalar are given, respectively, by

0 = p θ -ǫ θ , (6)
and

ω kθ sin(ωt + φ θ ) = ǫθ cos(ωt + φ ǫ θ ). (7)
An interesting feature of this last equation is that it does not contain any source term. On the level of single point averages it is therefore difficult to do any predictions on the behavior of kθ and ǫθ . The only result obtained from this equation is that

ω| kθ | = |ǫ θ |, (8)
and that the phase-shifts will differ by a constant value of π/2.

If we define the velocity and scalar timescales as

T = k ǫ , T θ = k θ ǫ θ , (9)
then we can relate the mixing rate to the scalar timescale by

χ θ = T -1 θ . (10)
Equivalently we can relate the energy transfer rate to the inverse of the integral timescale,

χ = T -1 . (11)
In the present case of statistically homogeneous flow, these definitions of the transfer and mixing rate are straightforward choices. In inhomogeneous flows the mixing process will also involve the turbulent diffusion and the definition of a mixing rate is less straightforward. When the scalar quantity is not passive, the definitions of mixing efficiency and mixing rate are far more complicated (as for instance in stably stratified flows [START_REF] Winters | Available potential energy and mixing in density-stratified fluids[END_REF], or flows generated by Rayleigh-Taylor instabilities [START_REF] Lawrie | Rayleigh-Taylor mixing in an otherwise stable stratification[END_REF] or Richtmyer-Meshkov instabilities [START_REF] Zhou | Asymptotic behavior of the mixed mass in Rayleigh-Taylor and Richtmyer-Meshkov instability induced flows[END_REF]).

Spectral description of the problem

The case we consider consists in an isotropic turbulent flow, maintained statistically stationary by a forcing containing a time-periodic component.

More precisely, we consider the dynamics of a turbulent flow, characterized by the kinetic energy spectrum E(κ, t), given by the evolution equation

∂ t E = -∂ κ Π -2νκ 2 E + P. (12)
where

P = P + P cos(ωt), ( 13
)
and ν is the kinematic viscosity. The choice of a model for the spectral flux Π determines the precise dynamics of the flow. An important property of this flux is that it should vanish at k = 0 and k = ∞, so that the integral of ( 12) yields (1).

We consider the advection of a passive scalar in the same flow, which is forced by a large scale injection, without any periodic component (P θ = P θ ).

The scalar spectrum E θ (κ, t) obeys,

∂ t E θ = -∂ κ Π θ -2ακ 2 E θ + P θ , ( 14
)
where α is the diffusivity of the scalar. The case of unity Prandtl number,

ν = α is investigated.
We consider that both the velocity and scalar forcing are confined to the largest scales. To facilitate the considerations, we assume the forcing terms to act on wavenumber κ f only,

P = p δ(κ -κ f ), P θ = p θ δ(κ -κ f ). (15)
A more complicated case could be investigated, where the periodic forcing acts at another scale than the steady part of the forcing, with possible relevance to flow-control, or geophysical flows, but this case is not considered here.

The energy spectrum and scalar spectrum will contain a time-averaged component and a periodic component,

E = E + Ẽ cos(ωt + Φ E ), (16)
E θ = E θ + Ẽθ cos(ωt + Φ E θ ). ( 17
)
The time-averaged spectra obey then the equations,

0 = -∂ κ Π -2νκ 2 E + P (18) 0 = -∂ κ Π θ -2ακ 2 E θ + P θ , (19)
and the periodic contributions,

∂ t + 2νκ 2 Ẽ cos(ωt + Φ E ) = -∂ κ Π cos(ωt + Φ Π ) + P cos(ωt)(20) ∂ t + 2ακ 2 Ẽθ cos(ωt + Φ E θ ) = -∂ κ Πθ cos(ωt + Φ Π θ ) . (21)
Upto this point we have not introduced any serious approximations, apart from the assumption that the response of the fields will be periodic with frequency ω. In the following, section 3 we will first focus on the low frequency dynamics of the velocity and scalar statistics. Then, in section 4, the high frequency dynamics of the scalar will be investigated.

Low frequency dynamics

In the low frequency, or quasi-static limit all scales have sufficient time to adapt to the forcing. The modulation does therefore not introduce any phase-shifts in the statistics. This means that in equations ( 20) and (21) the quantities Φ E , Φ E θ , Φ Π and Φ Π θ are close to zero. Also, the time-derivative vanishes in this limit. The equations simplify then to

2νκ 2 Ẽ = -∂ κ Π + P (22) 2ακ 2 Ẽθ = -∂ κ Πθ . (23)
To solve these equations, we need to specify the fluxes Π and Π θ . The simplest model for those fluxes, compatible with conserved flux and Kolmogorov-Obukhov scaling, is Kovaznay's closure for the nonlinear transfer of kinetic energy [START_REF] Kovaznay | Spectrum of locally isotropic turbulence[END_REF] and its analog for scalar turbulence [START_REF] Rubinstein | Reassessment of the classical closures for scalar turbulence[END_REF]. These models read

Π = C -3/2 k κ 5/2 E 3/2 , ( 24
)
with C k the Kolmogorov constant, and

Π θ = C θ κ 5/2 E 1/2 E θ , (25)
where C θ is related to the Kolmogorov and Corrsin-Obukhov constants by the relation

C θ = (C k C 2 CO ) -1/2 . ( 26
)
These closures are consistent with the dominant effect of straining in the transfer of energy and scalar variance towards small scales, compatible with Kolmogorov and Corrsin-Obukhov scaling. The influence of sweeping is not dominant in the energy transfer process, but its effect on the mixing efficiency constitutes an interesting direction for further research [START_REF] Zhou | Colloquium: Magnetohydrodynamic turbulence and time scales in astrophysical and space plasmas[END_REF][START_REF] Zhou | Renormalization group theory for fluid and plasma turbulence[END_REF]. Now these models are specified we can solve the equations to assess the influence of the modulated forcing on the dynamics of E and E θ .

Modulated kinetic energy

Introducing E = E + Ẽ in (24) and expanding in terms of Ẽ/ Ē we obtain,

Π ≈ C -3/2 k κ 5/2 E 3/2 1 + 3 2 Ẽ cos(ωt) E + O(( Ẽ/E) 2 ) , (27)
So that

Π ≈ C -3/2 k κ 5/2 E 3/2 1 + O(( Ẽ/E) 2 ) , (28)
and

Π cos(ωt) ≈ C -3/2 k κ 5/2 E 3/2 3 2 Ẽ cos(ωt) E + O(( Ẽ/E) 2 ) . (29)
Combining these expressions, we have therefore to first order

Π = 3 2 Ẽ E Π. ( 30
)
In the inertial range, in the quasi-static limit, fluxes should balance the energy input,

Π ≈ k P dk = p, Π ≈ k P dk = p, (31)
which gives

Ẽ = 2 3 p p E. (32)
Integration gives the modulated kinetic energy in the quasi-static limit,

k = 2 3 p p k, (33)
a result which was also obtained in [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] by elementary arguments.

Average kinetic energy and integral timescale

From equation (28) it can be seen that, to first order, the modulated energy input does not influence the time-averaged dynamics. This is consistent with previous studies of modulated turbulence, where the amplitude of the modulated forcing was small compared to the averaged component [START_REF] Kuczaj | Response maxima in timemodulated turbulence: Direct numerical simulations[END_REF], [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF].

However, in the recent simulations by Yang et al. [START_REF] Yang | Mixing in modulated turbulence: numerical results[END_REF], where the amplitude of the modulation was taken equal to the value of the time-averaged component, a significant modification of the averaged kinetic energy was observed for low frequencies. For such amplitudes, it is possible that the linear approximation breaks down, and higher order contributions should be taken into account.

Doing so, to second order, expression (28) becomes

Π = C -3/2 k κ 5/2 E 3/2 1 + 3 8 Ẽ2 cos 2 (ωt) E 2 . ( 34
)
The time-average of the quadratic fluctuation is,

Ẽ2 cos 2 (ωt) = 1 2 Ẽ2 , (35)
and combining this with (32), we obtain in the inertial range

p = ǫ = Π = C -3/2 k κ 5/2 E 3/2 1 + 1 12 p p 2 . ( 36
)
This shows that the modulation leads to an enhanced energy transfer to small scales. The physical consequence of this is that if there is a time-periodic part in the large-scale forcing of a flow, the energy spectrum will be equal to

E(k) = C ′ k ǫ 2/3 k -5/3 , ( 37
)
where the Kolmogorov constant is affected by the second-order correction.

Its effective value will be

C ′ k = C k 1 + 1 12 p p 2 -2/3 , (38)
and for a fixed energy input, the kinetic energy will be lowered by the same factor, with respect to the unperturbed value

k ∞ k = k ∞ 1 + 1 12 p p 2 -2/3
.

(39)

Modulated scalar variance

Introducing E = E + Ẽ and E θ = E θ + Ẽθ in (25) and expanding in terms of Ẽ/ Ē and Ẽθ / Ēθ we obtain to first order

Πθ ≈ C θ κ 5/2 E 1/2 E θ Ẽθ E θ + 1 2 Ẽ E . ( 40
)
Since there is no modulated production term, the flux in the inertial range should vanish. This gives the relation Ẽθ

E θ + 1 2 Ẽ E = 0, (41)
and thereby

Ẽθ = - 1 2 Ẽ E E θ . (42)
Using (32), we find that Ẽθ = -

1 3 p p E θ . ( 43
)
The minus sign in this expression shows that the scalar modulation is always completely out of phase with respect to the kinetic energy. The modulated scalar variance is obtained by integration, yielding

| kθ | = 1 3 p p k θ , (44)
and according to [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF],

|ǫ θ | = 1 3 p p ǫ θ ω. (45)

Average scalar variance and mixing rate

Again, as for the average kinetic energy, the modulation of the flow does not affect the linear expansion of the average scalar flux Π θ . Retaining second order contributions in the expansion of (25), we obtain

Π θ = C θ κ 5/2 E θ E 1/2 1 + 1 2 Ẽθ Ẽ E θ E cos 2 (ωt) - 1 8 Ẽ2 E 2 cos 2 (ωt) (46)
using ( 43), ( 35) and (32) this gives

Π θ = C θ κ 5/2 E θ E 1/2 1 - 1 12 p p 2 . (47)
Again, as for the kinetic energy, we see that the modulation affects the flux, but this time the effect is the opposite: the transfer rate is diminished. For a fixed input, in the inertial range we have,

p θ = ǫ θ = Π θ , (48)
yielding, with (37) and 38,

E θ = C ′ CO ǫ θ ǫ -1/3 k -5/3 , (49)
where

C ′ CO = C CO 1 + 1 12 p p 2 1/3 1 -1 12 p p 2 . ( 50
)
The variance of the passive scalar is then enhanced by the same factor,

k θ = k θ (p = 0) C ′ CO C CO . (51)
Since the mixing rate can be defined as

χ θ = ǫ θ /k θ , (52)
we have for a fixed scalar input,

χ θ (p) = χ θ (p = 0) C CO C ′ CO . ( 53
)
This expression is shown in Figure 1, together with the kinetic energy transfer rate. The mixing rate goes down for a large-amplitude modulation and the energy transfer rate increases, as observed in the DNS [START_REF] Yang | Mixing in modulated turbulence: numerical results[END_REF].

High frequency behavior of the modulated variances

In the previous section, the quasi-static limit was considered. In that limit the different phase-shifts disappeared from the dynamics. When the frequency goes up, the phase-shifts Φ E and Φ E θ , both functions of ω and κ, complicate the analysis. In [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF], it was shown that it is convenient to rewrite these equations in terms of spectra in-and out-of-phase,

E = E + F cos(ωt) + G sin(ωt), (54)
E θ = E θ + Fθ cos(ωt) + Gθ sin(ωt). ( 55
)
We have hereby eliminated the phases Φ E and Φ E θ from our description.

Both descriptions are obviously equivalent and can be related by

F 2 + G2 = Ẽ2 and tan(Φ E ) = -G/ F , (56)
with equivalent expressions for the scalar case. The first question to address is how the quantities Fθ and Gθ will be influenced by the periodically fluctuating energy input. The expressions for F and G were determined in [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF]. It was shown that in the inertial range, for large ω, F , G were given by

G/E ∼ ω -3 κ 2 f p and F /E ∼ ω -2 κ 2 f pǫ -1/3 κ -2/3 . (57)
In the light of our recent investigation of out-of-equilibrium turbulence [START_REF] Bos | Dissipation in unsteady turbulence[END_REF],

G corresponds to the equilibrium and F to the non-equilibrium contribution of the modulated kinetic energy spectrum.

The derivation of the large-frequency behaviour of kθ is straightforward.

First we replace in Π θ (expression (25)) the energy spectrum and the scalar spectrum by their expressions (54) and (55). Assuming the fluctuating spectra small compared to the time-averaged spectra, Ẽ/E ≪ 1, Ẽθ /E θ ≪ 1, the flux can be approximated (to first order) by, Πθ cos(ωt + φ Π ) ≈ Π θ Fθ cos(ωt) + Gθ sin(ωt)

E θ + 1 2 F cos(ωt) + G sin(ωt) E . (58)
We substitute this expression in the equation for the fluctuating scalar spectrum (21). Evaluating the expression at t = 0 and t = π/2 yields two coupled equations,

-ω Fθ = -∂ κ Π θ Gθ E θ + 1 2 G E -2ακ 2 Gθ (59) ω Gθ = -∂ κ Π θ Fθ E θ + 1 2 F E -2ακ 2 Fθ . (60)
In the inertial range, we assume that

E(k) = C k ǫ 2/3 κ -5/3 , E θ = C CO ǫ θ ǫ -1/3 κ -5/3 . (61)
We further use the results (57), and, ignoring the diffusive terms, the above set of equations can be rewritten as Fθ

E θ = ψ∂ κ Gθ E θ (62)
and

[1 + ψ∂ κ [ψ∂ κ ]] Gθ E θ = ψω -2 κ 2 f pǫ -1/3 κ -5/3 (63)
where we introduced the expression,

ψ ∼ ω -1 ǫ 1/3 k 5/3 . (64)
The second term in (63), containing the second derivative, is proportional to ω -2 . For large frequencies we can thus ignore it compared to the first one, yielding, Gθ

E θ ∼ ω -3 κ 2 f p (65) and Fθ ∼ ǫ 1/3 k 2/3 ω Gθ . (66)
Since p = ǭ ∼ k/T and κ f ∼ ǭk -3/2 , this can be expressed as Gθ

E θ ∼ (ωT ) -3 (p/p), (67)
and therefore kθ ∼ kθ (ωT ) -3 (p/p).

We see thus that the periodic part of the scalar variance, at large ω is determined by large scale quantities. Furthermore, the scalar dissipation is given directly by

|ǫ θ | = ω| kθ |, (69)
and therefore,

ǫθ = p p ǫ θ (ωT ) -2 . ( 70
)
Our estimates for the low and high frequency asymptotes are sketched in

Discussion

The main result obtained in the present manuscript involves the mixing rate. According to Figure 1, the mixing rate χ θ will go down to 89% of its value if the flow-forcing contains a modulation with an amplitude ratio p/p = 1. The energy-transfer rate, on the other hand will increase towards χ/χ(p = 0) = 1.05. These figures are the same order of magnitude as the figures obtained in the numerical study [START_REF] Yang | Mixing in modulated turbulence: numerical results[END_REF]. Note that at an amplitude of p p = 0.2 the change in the mixing rate and energy transfer becomes less than 1%. It is thus clearly a large amplitude effect. For smaller forcing amplitudes the linear prediction that the average fields are unaffected is approximately satisfied.

The same conclusions hold for the Kolmogorov and Corrsin-Obukhov constant. Their variations also become important only for large amplitude ratios. A recent related study of a similar effect can be found in [START_REF] Chien | Effects of fluctuating energy input on the small scales in turbulence[END_REF]. They exploited experimentally the idea by Monin and Yaglom [START_REF] Monin | Statistical Fluid Mechanics II[END_REF], section 25.1, who introduced a model to illustrate the influence of a fluctuating dissipation rate, or energy injection rate, on the Kolmogorov constant.

We can summarize the new physical insights as follows. The time-averaged nonlinear transfer is, at second order, influenced by the modulation of the energy injection. This enhances the energy flux. For a given energy injection at a given scale, this will decrease the level of the kinetic energy, or equivalently, the value of the Kolmogorov constant. The transfer of the passive scalar is also affected by the modulation, but the influence is the contrary. Qualitatively, this can be understood, at least partially, by the fact that the kinetic energy is decreased, so that the mixing of the scalar by the velocity fluctuations is less efficient. Adding a modulation to the forcing of a turbulent flow will thus decrease the integral velocity timescale and, simultaneously, increase the scalar timescale.

Figure 1 :

 1 Figure 1: Transfer rate χ and mixing rate χ θ in the low frequency limit, as a function of the relative forcing amplitude. Adding a slow, large amplitude modulation to the energy input diminishes the mixing rate and enhances the kinetic energy transfer rate.

Figure 2 :

 2 Figure 2: Frequency response of the modulated part of the passive scalar variance and dissipation.