Marc Budinger

Jean-Charles Passieux

Christian Gogu

Amine Fraj

Keywords: metamodels, surrogate modelling, scaling laws, Buckingham theorem, components sizing, preliminary/conceptual design

Scaling-law-based metamodels for the sizing of mechatronic systems

INTRODUCTION

A mechatronic system [START_REF]Design methodology for mechatronic systems[END_REF] [2] expands the capabilities of conventional mechanical systems through the integration of different technological areas (Figure 1) around:

• a power transmission part, which is a combination of components from mechanical, electrical, power electronics or fluid power technologies; • an information processing part, which is a combination of electronics, instrumentation, automatic, signal processing, and information technologies. 1 Corresponding author: Tel.: (+33) 561559960; Fax: (+33) 561559950; E-mail: marc.budinger@insa-toulouse.fr

The references [START_REF] Van Der Auweraer | Virtual engineering at work: the challenges for designing mechatronic products[END_REF] and [START_REF] Hehenberger | Hierarchical design models in the mechatronic product development process of synchronous machines[END_REF] highlight that the design of such multi-domain systems require different modeling layers as represented Figure 2: • A mechatronic layer, to take into account the functional and physical coupling between components. This level of modeling is usually done using 0D-1D models [START_REF] Van Der Auweraer | Virtual engineering at work: the challenges for designing mechatronic products[END_REF] also called lumped parameter models represented by algebraic equations, ordinary differential equations (ODE) or differential algebraic equations (DAE) [START_REF] Cellier | Continuous System Modeling[END_REF]. • A specific domain layer, to describe the performance limits and parameters necessary in the previous layer, based on a geometric representation. The specific domain phenomena are generally represented through partial differential equations (PDE). This level of modeling can be achieved, for simplified geometries using analytical models or, for complex 2D and 3D geometries, using numerical approximations like finite element method (FEM) for instance.

Figure 2 -Hierarchical design models (based on reference [4])

The design of the power part with a system integrator's point of view should allow to optimally size and specify the components of multiple technologies [START_REF] Roos | Optimal selection of motor and gearhead in mechatronic applications[END_REF] interacting together. This system level design, distinct of component design, needs to represent in the mechatronic layer the key informations of the specific domain layer with dedicated models [START_REF] Raymer | Aircraft Design: A Conceptual Approach[END_REF] [8] [START_REF] Liscouet | Modelling approach for the simulation-based preliminary design of power transmissions[END_REF]. The latter enables the designer to take multiple design constraints into account easily. They are referred to as "estimation models" in this paper. They directly and explicitly link a few primary characteristics, such as overall dimensions of components, to the secondary characteristics needed for the sizing [START_REF] Liscouet | Modelling approach for the simulation-based preliminary design of power transmissions[END_REF] and optimization [START_REF] Budinger | Optimal preliminary design of electromechanical actuators[END_REF]. The capacities required of these estimation models are: to present a form that is simple to handle and to implement in different calculation tools, to lend themselves to easy analytical manipulations, and to be reusable in an area slightly different from the one where they were initially employed. To satisfy these constraints, simplified analytical models are often used [START_REF] Gomis-Bellmunt | Hydraulic actuator modeling for optimization of mechatronic and adaptronic systems[END_REF] [START_REF] Gomis-Bellmunt | Linear electromagnetic actuator modeling for optimization of mechatronic and adaptronic systems[END_REF]. Among them, scaling laws have proved effective to represent a physical phenomenon over wide ranges of variation [START_REF] Pahl | Engineering Design: A Systematic Approach[END_REF] [START_REF] Budinger | Estimation models for the preliminary design of electromechanical actuators[END_REF]. However, these models are valid under certain conditions, among which one can mention geometry and material similarities, and uniqueness of the driving physical phenomenon. For the system designer the models should be as predictive as possible. Detailed finite element models, able to precisely predict the physical phenomena are still too time-consuming in such a context. Despite a recent thrust of work on model order reduction [START_REF] Bialecki | Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis[END_REF] [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF], the computational cost of finite element models remains prohibitive in the preliminary design phase. The use of metamodelling techniques [START_REF] Simpson | Metamodels for Computer-Based Engineering Design: Survey and Recommendations[END_REF] [17] [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF] is thus interesting for this purpose. In this paper, we propose an enhanced power law type model based on metamodels to represent the models of the specific domain layer (Figure 2) into the mechatronic layer.

The second section of this paper will show the interest of scaling laws to establish the estimation models needed by the system designer. In the third section, a metamodeling method based on scaling laws will be proposed to extract simple, global expressions of estimation models from local numerical simulations (FEM). The regression process is first illustrated with a mathematical function and then with two examples of mechatronic components from different domains: a limited-angle electromagnetic actuator and a flexible mechanical hinge. Noise effects and comparisons with classical metamodelling techniques are illustrated through these examples.

ESTIMATION MODELS WITH SCALING LAWS

Scaling laws and Buckingham theorem

Scaling laws based on dimensional analysis, also called similarity laws or allometric models, have been very successfully used throughout the past decades for solving scientific and engineering problems and for presenting results in a compact form. In the design of mechatronic systems scaling laws [START_REF] Budinger | Estimation models for the preliminary design of electromechanical actuators[END_REF] [START_REF] Pahl | Engineering Design: A Systematic Approach[END_REF], allow estimation models to be obtained from a single reference component by using three main modeling assumptions: a. Material similarity: all material and physical properties are assumed to be identical to those of the component as the reference;

b. Geometric similarity: the ratio of all the lengths of the component under consideration to all the lengths of the reference component is constant;

c. Uniqueness of design driver: only one main dominant physical phenomenon drives the evolution of the secondary characteristic y.

The mathematical form of scaling laws is a power law: This form and the conditions for obtaining equation ( 1) are demonstrated here using the Buckingham theorem [START_REF] Taylor | Dimensional Analysis for Engineers[END_REF] [20] [START_REF] Miragliotta | The power of dimensional analysis in production systems design[END_REF]. An estimation model seeks to identify a relationship between 2+n+m parameters:

ݕ = ܮ݇ (
f ቆy, L, d ଵ , d ଶ , … , d ୬ ᇣᇧᇧᇧᇤᇧᇧᇧᇥ ଵା୬ , p ଵ , … , p ୫ ᇣᇧ ᇧᇤᇧ ᇧᇥ ୫ ቇ = 0 (2)
with:

• 1 parameter corresponding to the secondary characteristic y to be estimated;

• 1+n parameters characterizing the geometrical dimensions L and d i ;

• m parameters characterizing physical and material properties p i ;

Depending on the number of physical units u (e.g. m, kg, s, etc.) involved in the problem, this relationship can be rewritten using dimensionless parameters ߨ :

݂′ ൭ߨ ௬ , ߨ ଵ , ߨ ଶ , … , ߨ ᇣᇧ ᇧ ᇧᇤᇧ ᇧ ᇧᇥ , ߨ ଵ , … , ߨ ᇱ ᇣᇧ ᇧ ᇧᇤᇧ ᇧ ᇧᇥ ᇲ ൱ = 0 (3)
where:

ߨ ௬ = ܮݕ ಽ ෑ  (4) ߨ = ݀ ܮ (5)
ߨ = ܮ ಽ ෑ  ೕ ( 6
)
The number m', usually smaller than m, depends on the number of physical units u as expressed by the Buckingham theorem:

݉ ᇱ = 1 + ݉ -ݑ (7
) If only one main simple physical phenomenon drives the evolution of the secondary characteristic y, the number m' is often equal to zero. If m' is not equal to zero, the remaining dimensionless numbers can generally be expressed through ratios of material properties with similar units. These m' dimensionless number thus do not depend on L. With this condition and the first two assumptions a. and b., the π i and π pi dimensionless numbers are constant and it follows that:

ߨ ௬ = ܮݕ ಽ ෑ  = constant (8)
which gives relation (1) if we assume the material or physical properties p i to be constant during the sizing.

Examples of scaling laws

To illustrate the construction and use of such laws, two examples will be given here. These examples address conventional components of mechatronic systems: brushless motors and bearings. We assume here that the main design criterion for the motor is the winding temperature. The dominant thermal phenomenon will be assumed to be convective. The bearings are for their part designed to withstand a maximum mechanical stress.

For the brushless motors, equations (2), ( 3) and (1) become the following:

• For thermal aspects: the current density J can be linked to the dimensions through

݂ ൭,ܬ ,ܮ ݀ ଵ , ݀ ଶ , … , ݀ ᇣᇧᇧᇧᇤᇧᇧᇧᇥ ଵା , ߩ, ߠ, ℎ ᇣᇤᇥ ൱ = 0 (9)
which, according to the Buckingham theorem, leads to

݂′ ቌ ఘ మ ఏ , ௗ భ , ௗ మ , … , ௗ ᇣᇧᇧᇤᇧᇧᇥ ቍ = 0 where m'=0 (10)
and

ܬ = ܮ݇ ି,ହ (11)
• For magnetic aspects : the torque T can be linked to dimensions and current density through

݂ ൭ܶ, ,ܮ ݀ ଵ , ݀ ଶ , … , ݀ ᇣᇧᇧᇧᇤᇧᇧᇧᇥ ଵା , ,ܬ ܤ ถ ൱ = 0 (12)
which, according to the Buckingham theorem, leads to

݂′ ቌ ் ర , ௗ భ , ௗ మ , … , ௗ ᇣᇧᇧᇤᇧᇧᇥ ቍ = 0 where m'=0 (13)
and

ܶ = ܮܬ݇ ସ (14)
with : J the current density, L the length of the motor, d i other geometrical dimensions, ρ the resistivity of the copper, θ the maximal temperature rise for the winding insulation, h the convection coefficient, B r the remanent induction of the permanent magnet, and T the electromagnetic torque.

Combining these two aspects provides an estimate of torque depending on motor size:

ܶ = ܮ݇ ଷ,ହ (15)
The assumptions a. and b. allow the motor weight to be estimated:

ܯ = ܮ݇ ଷ = ݇ܶ ଷ/ଷ,ହ (16)
With the same approach, the weight of the bearing can be estimated from their load-bearing capacity C: ܥ = ܮ݇ ଶ and ܯ = ܮ݇ ଷ = ܥ݇ ଷ/ଶ (17) More details and examples of the development of such laws for mechanical and electromechanical components are given in [START_REF] Budinger | Estimation models for the preliminary design of electromechanical actuators[END_REF]. Figure 3 and Figure 4 compare these relationships to data from industrial catalogues, and show that scaling laws can provide good fits for the quantities of interest of such components

Interests and limitations of scaling laws

Scaling laws have assets that make them attractive for the design of mechatronic systems [START_REF] Pfenning | Recent Advances Towards an Integrated and Optimized Design of High-lift Actuation Systems[END_REF] [23] [START_REF] Liscouet | Modelling approach for the simulation-based preliminary design of power transmissions[END_REF]. Their simple form makes them easy to manipulate and customize as they require only one reference to determine the multiplier coefficient k, the power coefficient a being determined by the physical phenomena. They have a monotonous progression valid over a wide range of sizes (several orders of magnitude) which avoids the risk of possible mathematical aberrations of metamodels used outside their construction bounds.

However they have some limitations. Although the similarity of the materials can be easily verified for a given technology, the geometric similarity is not necessarily verified or sought. For motors (Figure 3), this point mainly explains the estimation errors of scaling laws in Figure 3. Manufacturers have a tendency to use the same motor diameter (i.e. the same iron sheet) for different lengths and torques. Obtaining a scaling law also requires a dominant physical phenomenon. For example, in the case of the electric motor, if the conductive heat transfer phenomenon cannot be neglected, relations [START_REF] Liscouet | Modelling approach for the simulation-based preliminary design of power transmissions[END_REF] and [START_REF] Budinger | Optimal preliminary design of electromechanical actuators[END_REF] should take the form:

݂ ൭,ܬ ,ܮ ݀ ଵ , ݀ ଶ , … , ݀ ᇣᇧᇧᇧᇤᇧᇧᇧᇥ ଵା , ߩ, ߠ, ℎ, ߣ ଵ , ߣ ଶ , … ᇣᇧᇧ ᇧᇧᇤᇧᇧᇧᇧᇥ ൱ = 0 (18)
which gives

݂′ ቌ ఘ మ ఏ , ௗ భ , ௗ మ , … , ௗ ᇣᇧᇧᇤᇧᇧᇥ , ఒ భ , ఒ మ , … ᇣᇧᇤᇧᇥ ᇱ ቍ = 0 (19)
with λ i the thermal conductivities of the various materials of the motor.

In this case, either a global scaling law would have been impossible to establish or its domain of validity would have been smaller.

SCALING LAWS BASED ON REGRESSION

The objective of this section is to show how it is possible to increase the use and validity of scaling laws through appropriate regression techniques as proposed by our method.

Form of models

To maintain the physical meaning and benefits of scaling laws, the regression model will be based on the power form [START_REF]Design methodology for mechatronic systems[END_REF]. However, to eliminate assumptions b. (geometric similarity) and c. (uniqueness of design driver), estimation models of the following, more general form will be sought:

ݕ = ,ܮ(݂ ߨ ଵ , ߨ ଶ , … ) = ݇(ߨ ଵ , ߨ ଶ , … ܮ) (గ భ ,గ మ ,… ) with ߨ = ௗ (20)
where :

• y is the parameter to be estimated;

• L is the main dimension of the system;

• d i are the secondary dimensions of the system;

• π i are the dimensionless numbers representative of form factors of the system.

With geometrical similarity, the form factors are constant (π i =C st ) and the form (20) simplifies into the classic scaling law form [START_REF]Design methodology for mechatronic systems[END_REF]. The desired form of model ( 20) lends itself to direct regression less easily than polynomial response surfaces [START_REF] Raymond | Response surface methodology : process and product optimization using designed experiments[END_REF], radial basis functions [START_REF] Simpson | Metamodels for Computer-Based Engineering Design: Survey and Recommendations[END_REF] [25], or kriging [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] [27] directly on the parameter y. The objective of the coming sections is to give an approach to determine the shapes of functions ݇(ߨ ଵ , ߨ ଶ , … ) and ܽ(ߨ ଵ , ߨ ଶ , … ) representative of multiplier and power coefficients.

Regression process

Figure 3 describes the overall meta-modelling process proposed to carry out regressions of form [START_REF] Murphy | Similitude in engineering[END_REF]. The three main steps of the process are described below and the details of each step are described in the following subsections.

• Step 1 -data generation: the objective is first to generate the combinations of variables L and π i to be simulated. With deterministic computer experiments, these sample points should be chosen to fill the design space. Once the DoE is defined, a call is then made to the finite element codes, or any other sizing code, to obtain the desired data y. The data y are then projected according to form (20) by calculating the coefficients k and a for each configuration π i of form factor parameters.

• Step 2 -study of the general shape of the function: the evolution of each coefficient k and a is then analysed to determine the overall shape to be imposed on function [START_REF] Murphy | Similitude in engineering[END_REF] through the choice of the forms of functions ݇(ߨ ଵ , ߨ ଶ , … ) and ܽ(ߨ ଵ , ߨ ଶ , … ). These forms of functions cannot be initially assumed and depend on each problem. • Step 3 -building of the metamodel: based on the function forms selected in the previous step, a regression process is applied to simulation data. A comparison of the results predicted using Eq. ( 20) with the initial data y can then validate the quality of these regressions.

Even though the model form of Eq. ( 20) cannot be expressed by linear combination of basis functions, some of its characteristics are interesting and will be used for the proposed meta-modeling process. The following subsections describe the choices made here for each of these steps according to these characteristics.

Step 1 -Data generation

The range of variation for the parameter L is potentially very large and is potentially much smaller for the shape parameters π i . The proposed process will use a DoE (design of experiments) which is logarithmically distributed on n L levels for the parameter L and linearly distributed on n π levels for the p parameters π i . This simple DoE enables the design space to be filled. This DoE (step 1.1), mainly similar to a full factorial design, requires at least three levels for each parameter in order to calculate (steps 1.3) and analyse (steps 2.1 and 2.2) coefficients k and a.

According to this DOE, calls are then made (step 1.2) to the finite element codes, or any other sizing code, to obtain the desired data y. The final number of data y is n L . n π p . The conventional regression techniques use these simulation data directly. Here in order to adapt the data to form [START_REF] Murphy | Similitude in engineering[END_REF] and as a power law can be represented linearly in logarithmic-scale graphs, the projection of the data y onto the coefficients a and k (step 1.3 of Figure 5) is done in a logarithmic plane plotting the y according to L data for each combinations of the form factors π i . Figure 6 represents the graphs where, for each combination of coefficients π i ., the evolution of y according to L follows a power law and is therefore described by a straight line in a logarithmic plane. The algorithm computes the slopes a of the segments between two consecutive values of L and returns the mean average slope mean(a) for each combinations of π i . The multiplier coefficient k is calculated for each curve by minimizing the least square error for the average slope. In order to check that the data are well fitted by power laws, a quality criterion is also calculated. This criterion expressed by the relation std(a)/mean(a) is usually very close to zero and shows that the problem is well expressed by using a scaling law for a given configuration of form factors. After this step, the n L .n π p values y are transformed into n π p values vectors for all coefficients a and k.

.2. Step 2 -Metamodel form definition

The objective of this step is to determine the form that can be given to functions ܽ(ߨ ଵ , ߨ ଶ , … ) and ݇(ߨ ଵ , ߨ ଶ , … ).

The power coefficient ܽ(ߨ ଵ , ߨ ଶ , … ) of scaling laws is constant if the physical phenomenon that drives the design of the system does not depend on the form factors π i . If the physical phenomenon switches from one to another (e.g. convection to conduction for thermal transfer), the power coefficient is a function that changes between 2 extreme values. If the switch between the two physical phenomena is progressive (often the case for components of mechatronics systems) then the surface ܽ(ߨ ଵ , ߨ ଶ , … ) will be characterized by a small number of peaks and valleys and can potentially be well approximated by a polynomial function. It remains to determine the exact form of the polynomial function. Since all the parameters π i do not usually affect the power coefficient it is useful to conduct a sensitivity analysis at this point. Step 2.2 thus corresponds to a sensitivity analysis, which is conducted here by performing a linear regression after normalization between -1 and +1 of parameters π i . The analysis of the regression coefficients leads to the quantification of the average influence of each parameter. Only relevant parameters will be considered in the final expression of ܽ(ߨ ଵ , ߨ ଶ , … ). Note that any other technique of sensitivity analysis could be used here (e.g. ANOVA [START_REF] Damon | Experimental design, ANOVA, and regression[END_REF] or Sobol's indices [START_REF] Sobol | Sensitivity estimates for non-linear mathematical models[END_REF]).

The multiplication coefficient ݇(ߨ ଵ , ߨ ଶ , … ) of scaling laws can often take the form of products of functions, some of which may also be expressed in the form of power laws. Step 2.1 (Figure 5) analyses these options by using a logarithmic plane again as the functions that can be decomposed into products are represented by parallel curves in logarithmic-scale graphs (Figure 7a). The power laws are represented by straight lines (Figure 7b). If:

݇(ߨ ଵ , ߨ ଶ , … ) = ݂(ߨ ଵ )݃(ߨ ଶ , … ) or ݈݃൫݇(ߨ ଵ , ߨ ଶ , … )൯ = ݈݃൫݂(ߨ ଵ )൯ + log(݃(ߨ ଶ , … )) (21)
the evolution of values of ݇(ߨ ଵ , ߨ ଶ , … ), a function of π 1 for different π i≠1 configurations, can give 3 possibilities: factorizable (Figure 7a), factorizable with a power law (Figure 7b), or non-factorizable (Figure 7c).

This type of analysis can be carried out for each variable π i and quantified by indicators describing: These criteria can be condensed graphically on a plane, as illustrated in Figure 7d, to quickly assess the form that the function ݇(ߨ ଵ , ߨ ଶ , … ) should have.

•

Step 3 -Metamodel building

The last step aims to build the metamodels from data generated in step 1 and according to the forms of functions defined in step 2. This process of regression is done in two sub-steps:

• The first uses the data a and k (outputs of step 1.3, Figure 5): a least squares polynomial regression, using the most influential parameters, is directly performed on the data a to obtain an approximation for ܽ(ߨ ଵ , ߨ ଶ , … ). Depending on the decomposition of ݇(ߨ ଵ , ߨ ଶ , … ) derived from the indicators of step 2.1 (Figure 7d), the regression is performed on each of the factorizable functions. Power laws are processed first by linear regression in logarithmic planes. Functions that cannot be expressed by power laws are approximated by polynomial functions using least squares regression. • The second uses y data directly (outputs of step 1.2, Figure 5): this nonlinear regression is initialized with results of the first regression and minimizes relative errors thanks to the lsqnonlin function of Matlab.

Finally a comparison of the two regression functions with respect to the initial data y can validate the quality of these regressions.

Example with a test function

A purely numerical example will help to illustrate and detail the various steps of the process described in Section 4.2. This test function, which uses the form of equation [START_REF] Murphy | Similitude in engineering[END_REF], is:

ݕ = ߨ ଵ ିଵ (4 + 5ߨ ଶ ߨ ଷ )ߨ ସ ଷ ܮ 2+0.1ߨ 2 (22)
A numerical uniformly distributed noise can be added to y data in order to test the robustness of the process to the possible computer noise of design codes or to variations due to neglected phenomena. For the first numerical application, this noise is set to +/-1% of y data. The different stages of regression, shown from Figure 8 to Figure 11,are: • Step 1.1, the generation of a design of experiments: L with 4 levels between 1 and 10 4 and each parameter π i with 3 levels between respectively .5 and 2, 1 and 10, 2 and 7, 1 and 10. The number of simulations and data y is 4.3 4 =324. • Step 1.3, the y data projection: The straight lines (cf. Figure 8) imply the power law nature of the relationship. For the example of Figure 8, the slope a is equal to 3. The deviation from a power law is expressed by a quality indicator, std(a)/mean(a) as explained in 3.2.1, which has the value of 0.24% and also appears in Figure 8. Without the parameter L the remaining design of experiments now has 3 4 =81 values for each parameter a and k.

•

Step 2, a study of the general shape of the function: Figure 9 shows the analysis (step 2.1) of the possible form for the function ݇(ߨ ଵ , ߨ ଶ , … ). From indicators presented in Figure 7 and section 4.2, it follows that the parameters π 1 and π 4 are suitable for factorization using power laws. Figure 10 illustrates the sensitivity analysis (step 2.2) to determine parameters influencing the evolution of the a coefficient. Only the parameter π 2 has to be kept here. We thus select an equation with the following form:

ݕ = ݇(ߨ ଵ , ߨ ଶ , … ܮ) (గ భ ,గ మ ,… ) = ݇′(ߨ ଶ , ߨ ଷ )ߨ ଵ భ ߨ ସ ర ܮ (గ మ ) (23)
with ݇′(ߨ ଶ , ߨ ଷ ) and ܽ(ߨ ଶ ) second order polynomial functions. • Step 3, the construction of the metamodel based on functions forms selected in the previous step: functions ݇′(ߨ ଶ , ߨ ଷ ) and ܽ(ߨ ଶ ) will be approximated by polynomial function. The first regression allows functions ݇(ߨ ଵ , ߨ ଶ , … ) and ܽ(ߨ ଵ , ߨ ଶ , … ) to be obtained independently by a log plane study (for the power law) or linear regression (for the other function). This first regression (step 3.1) gives : ݕ = (4.84 -0.14π ଶ -0.37π ଷ + 5.

01ߨ ଶ ߨ ଷ + 0.01π ଶ ଶ + 0.03π ଶ ଶ )ߨ ଵ ି.ଽଽ ߨ ସ ଷ ܮ (ଶା.ଵగ మ ) (24)
The second non-linear regression (step 3.2) working directly on y and π i data enables us to find : ݕ = (3.988 + 0.002π ଶ + 0.002π ଷ + 5.001ߨ ଶ ߨ ଷ )ߨ ଵ ିଵ ߨ ସ ଷ ܮ (ଶା.ଵగ మ ) (25) For the 2 regressions, the mean square error minimizes the relative and not the absolute errors. The regression quality check shows that less than 1% of error is achieved with equation ( 25). Without noise, exact relationship is found directly with the first and second regression. Even though numerical simulations are deterministic, errors may exist due to simplifications of the mathematical model or systematic errors of the numerical scheme. To quantify the effect of errors on the proposed process more precisely, we assume here that the error is similar to white noise of amplitude proportional to the quantity of interest y. The calculation process was repeated for different values of maximum noise: 0.2, 0.5, 1, 2 and 4%. Figure 9 shows the results of step 2.1 and 3.3 for these different levels. Figure 11a) was plotted with 10 runs for each noise level. Figure 11b) quantifies the maximum error found between the expression and the simulation results (with noise) or the pure analytical formula (without noise). From these results we can conclude that:

• Step 3 (regression), with its two regression sub-steps, is quite robust to noise and even decreases the influence of noises.

• Step 2.1 (determination of the shape of the coefficient k) is more sensitive to noise but allows the form of the function to be distinguished for noise levels less than + / -2%. The FEM simulations should therefore have a meshing quality sufficient to allow for this step.

Comparison with other analysis and regression methods

Reference [START_REF] Simpson | Metamodels for Computer-Based Engineering Design: Survey and Recommendations[END_REF] gives a classification of the different approximation techniques according to 3 criteria: the type of DoE, the choice of model and the model fitting approach. Here, the proposed meta-modelling approach is based on:

• a full factorial DoE with logarithmic distribution of levels distribution for the main dimension L parameter;

• a model choice based on scaling laws and power laws;

• a fitting approach using mainly linear and non-linear least square regression based on a projection of the data on scaling law coefficients.

Classic regression procedures generally start with a screening step mainly focusing on selection of the input parameters. The process assumes that this selection has been performed thanks to a previous sensitivity analysis or the designer's knowledge. Preliminary analysis, here, is devoted to determining the exact form of the model. Indicators and a graphical synthesis, similar to the Morris method [START_REF] Morris | Factorial Sampling Plans for Preliminary Computational Experiments[END_REF], have been proposed for this purpose.

Compared to the classical meta-modeling approaches (polynomial, kriging, radial basis functions, etc) the advantage of the method proposed here lies mainly in its physical justification, demonstrated here by using the Buckingham theorem. A recent review about dimensional analysis and pi-theorem illustrated by an example of application to the modeling of a Flexible Manufacturing System can be found in [START_REF] Miragliotta | The power of dimensional analysis in production systems design[END_REF]. Some other authors [START_REF] Gogu | Dimensionality Reduction Approach for Response Surface Approximations: Application to Thermal Design[END_REF] [32] [START_REF] Lacey | The Use of Dimensional Analysis to Augment Design of Experiments for Optimization and Robustification[END_REF] have also used this theorem with the main objective of reducing the number of parameters to be handled by polynomial regressions. Other authors [START_REF] Mendez | Scaling laws from statistical data and dimensional analysis[END_REF] have used dimensional analysis to determine constraints on power law regressions. Their approach, like other approaches solely based on dimensional analysis, however requires handling of all the parameters representing the dimensions and physical properties. Moreover, the scaling laws employed are pure power laws without the possibility of using other functions to express multipliers or power coefficients.

The approach proposed here does not directly use dimensional analysis but the particular form of the functions representing the physical characteristics of a component subjected to a change in size in the case of particular similarities. Previous studies [35] [20] on similarity have also worked on the development of models valid in the case of non-similarities (called distored models) but have not been applied to numerical simulations within a general framework. The possibility of keeping constant the ratios of some of the dimensions reduces the complexity of representation of the detailed design models. A compact model is then obtained with only the relevant design parameters for the design of mechatronic layer without having to handle all the parameters of the detailed design.

CASE STUDIES

The objective here is to illustrate the proposed regression process on mechatronic components from different domains. The main focus is placed on the illustration of the regression process and not to on the complete design study of a system. Reference [START_REF] Budinger | Optimal preliminary design of electromechanical actuators[END_REF] presents a methodology for the optimal preliminary design of electromechanical actuators using a model-based approach with different types of models (estimation, simulation, evaluation and meta-model).

The components considered are:

• an electromagnetic limited angle actuator [START_REF] Tsai | Design and control of a brushless DC limited-angle torque motor with its application to fuel control of small-scale gas turbine engines[END_REF], as shown in Figure 13a from [START_REF] Bencheikh | Limited angel torque actuator for fine angular positioning[END_REF], which illustrates the thermal and electromagnetic domains; • a flexural hinge composed of thin blades, as shown in Figure 13b from [START_REF] Santos | High accuracy flexural hinge development[END_REF], which illustrates the structural and mechanical domains.

The regression process could be applied to other types of motor, such as a brushless motor or a voice coil actuator [START_REF] Zhou | Design of Fast Steering Mirror systems for precision laser beams steering[END_REF], other flexural hinge [START_REF] Zhou | Design of Fast Steering Mirror systems for precision laser beams steering[END_REF] [START_REF] Santos | High accuracy flexural hinge development[END_REF] or any other physical mechatronic components. These 2 components associated with a mirror, a position sensor and control electronic could be used to build up a scanning mechanism [START_REF] Vinals | MTG scan mechanism predevelopment: design \& performance[END_REF]. The overall design of such mechatronic systems requires consideration of the various constraints that have a direct impact on the size of the components. These design constraints may come from different types of specifications:

• Performance specifications, particularly in terms of motion range, accuracy or bandwidth: the motion is linked to the rotational stiffness of the hinge and to the torque capability of the LAT actuator; the accuracy and the bandwidth are linked with resonance modes, which depend on the inertia and stiffness of components. • Endurance specifications, particularly in terms of life time and resistance to extreme stresses: the flexural link should be used at stress levels consistent with these specifications.

These two points are influenced by the environmental mechanical stress (vibration, acceleration) or thermal stress (a function of the type of heat transfer and the ambient temperature). Embedded applications add constraints of integration (size, volume or mass). All these constraints can be represented (Figure 12) by design scenarios such as: a) Thermal and mechanical mission profiles: these mission profiles are representative of the movement of the mirror during scanning. These time simulations can be used to calculate the average electromagnetic torque developed by the LAT. One can note that there exists a coupling of this torque with the rotational stiffness of flexural hinges.

b) Stability analysis: the transfer function of the overall multiphysic system can be used to estimate the closed loop performance especially in terms of bandwidth and stability margins. The indicators will be significantly influenced by the rotational resonance modes.

c) Vibratory environment: transient or frequential profiles for acceleration at the support can be used to represent the aggressive vibratory environment. The transverse or longitudinal modes of resonance are all influenced by the transverse or longitudinal stiffness of the flexible pivots and the mass of components.

Figure 12 -Estimation models : interface between mechatronic and domain specific layers

To manage all these constraints and iterate quickly between mechatronic and domain specific layers during preliminary sizing, it may be interesting to have estimation models expressing, for example:

• for the LAT actuator, torque or inertia depending on the dimensions.

• for the flexure hinge, stiffness (in translation or in rotation), maximal stress depending on the dimensions.

The following paragraphs show how these models can be obtained using the regression process presented in this article. The results of scaling-law-based metamodels are compared to polynomial approximations [START_REF] Raymond | Response surface methodology : process and product optimization using designed experiments[END_REF] or Radial Basis Functions (RBF) [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] and kriging [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF]. a) Limited Angle Torque (LAT) actuator (source : [START_REF] Bencheikh | Limited angel torque actuator for fine angular positioning[END_REF]): thermal and electromagnetic domains b) Flexural hinge (source [START_REF] Santos | High accuracy flexural hinge development[END_REF]): mechanical domain

Limited angle actuator study

The objective here was to illustrate how it is possible to link the continuous torque T to the main dimensions, the length, L, and the diameter, d, of the actuator. All radial dimensions were geometrically similar to d. With such a choice, one has two design parameters influencing significantly the torque and inertia of the actuator, two main characteristics for the system level design, without representing all dimensional parameters used during the component level design. The desired equation thus took the following form:

ܶ = ,ܮ(݂ ߨ ଵ ) = ݇(ߨ ଵ ܮ) (గ భ ) with ߨ ଵ = ௗ (26)
The design code, to which the regressions were applied, was set up using two 2D finite element simulations (Figure 14) with FEMM software [START_REF] Meeker | Finite Element Method Magnetics Software[END_REF]. A first thermal simulation found the permissible current density J for static continuous operation. The thermal environment was modeled by convective heat transfer between the periphery of the actuator and the environment. A second electromagnetic simulation used the previously calculated current density in order to obtain the continuous torque.

Several cases illustrate the regression method:

• The winding was moulded with or without a potting resin that ensured better conductive thermal transfer. This point addressed the effect of assumption c. (subsection 3.1). An equivalent thermal conductivity of copper and resin or air was deduced from a homogenization of the winding [START_REF] Idoughi | Thermal Model With Winding Homogenization and FIT Discretization for Stator Slot[END_REF] [44]. • Two DoE for different ranges were used to test the robustness and validity of the regressions obtained outside their original design. The two DoE had 4 levels for L and 6 levels for π 1 which led to 24 experiments. The difference between the DoE came from their range: for L, [20;60] Next, the proposed approach was applied to a winding with potting. The form of the scaling-law-based metamodel (SLAWMM) was:

ܶ = ݇(ߨ ଵ ܮ) (గ భ ) = ܿ ଵ ߨ ଵ ଶ ܮ య ା ర గ భ ା ఱ గ భ మ (27)
It was compared to three other approximations depending on the two design variables L and d: a 3 rd order polynomial response surface (RSM) [START_REF] Raymond | Response surface methodology : process and product optimization using designed experiments[END_REF], a Gaussian Radial Basis Function (RBF) [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] and a kriging approximation [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]. The polynomial functions were built with the same DoE as the scaling-law-based metamodels. The RBF and kriging approximations were built with a specific DoE: a latin hypercube (LHC) with the same range and number of experiments as the scaling-law-based metamodels, in this case 24. The polynomial regression, RBF and kriging functions were obtained with L and d parameters and not with L and π 1 which would not have been such a good choice. The least square regression for polynomial functions has minimized the relative error and not the absolute error in order to avoid gross errors on small actuators. The Gaussian width factor of RBF was chosen by minimizing a cross validation error [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]. The kriging approximation used a quadratic regression model and a Gaussian correlation function.

Figure 16 illustrates step 3.3, the model quality check, of the regression process. These validations were performed by comparing the 48 simulated results, 24 for the SLAWMM DoE and 24 for the LHC DoE, with the estimated results using each metamodel. Each approximation was thus tested with 24 configurations of the initial DoE that served for their construction and 24 other configurations within the same range. Note that the DoE using range 2 was much wider than when using range 1, leading to four additional orders of magnitude of variations on the output quantity that was being approximated. In all cases (range 1 and range 2), the scalinglaw-based metamodels (SLAWMM) showed good prediction qualities: less than 2 and 10% of maximum error for the range 1 (Figure 16a) and range 2 (Figure 16b) regression. The expressions obtained for equation [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF] differ slightly between the 2 ranges (see figure captions). Polynomial approximations, on the other hand, even with twice number of coefficients to be determined, showed a larger maximum error : 3% (Figure 16a) and 79% (Figure 16b) for range 1 and 2 respectively. For range 2, the maximum error was 18.5 % at the points of the DoE used for polynomial response surface construction but increases to 79% with the 24 configurations of the DoE that did not serve for its construction (predictive capability). Radial basis functions and kriging led to even worse results. These two metamodels had difficulty in coping with the very wide range of variation of the output quantity, which led, in particular to large errors for very small values. When the predictive capability of metamodels was assessed by testing them at points well outside their construction bounds, regression-based scaling laws remained valid (less than 15% error in Figure 16c). Gaussian RBF and kriging functions, with their high flexibility enabling a multiple peaks surface response to be represented, were at a disadvantage in this case of fairly monotonous surfaces: they could give very poor results at points involving large extrapolation. All three classic regression techniques (RSM, RBF and kriging) gave very poor results, even negative torque for small actuators, when they were used well outside the points that served for their construction (equations of range 1 used into range 2 as illustrated in Figure 16c). Power law forms are thus interesting as "structural components" for the specific metamodels of sizing problems where geometrical dimensions can have a wide range of variation. This robustness property can be important when conducting explorative design studies in areas not originally planned. As the regression form has a physical basis, some coefficients can be be fixed after a first complete study. For example, the c 3 coefficient of equation ( 27) can be fixed at 3.5 (see section 3.2) and future FEM simulations with different configurations have to find just 4 coefficients (c 1 , c 2 , c 4 and c 5 ) with a reduced size DOE.

ܶ = ݇(ߨ ଵ ܮ) (గ భ ) = 1.29.10 ି ߨ ଵ ଶ.ହଽ ܮ ଷ.ସ଼ି.ଶଵగ భ ା.ଷగ భ మ (T in N.m, L in mm) a)

Flexural hinge study

The objective here was to illustrate how it is possible to link the rotational stiffness, K r , to the main dimensions: L the length, d the diameter, e the thickness, and r the blade root radius (Figure 17a) of a flexural hinge composed of 2 thin blades of titanium alloy. The desired equation thus took the following form:

ܭ = ݂(݀, ߨ ଵ , ߨ ଶ , ߨ ଷ ) = ݇(ߨ ଵ , ߨ ଶ , ߨ ଷ )݀ (గ భ ,గ మ ,గ య ) with ߨ ଵ = ௗ , ߨ ଶ = ௗ and ߨ ଷ = ௗ (28)
The design code, to which the regressions were applied, is set up using a 3D finite element simulation (Figure 17b) with ANSYS software [45]. This static analysis gave the angular deflection for a given torque on the hinge, which, by ratio, gave the rotational stiffness. The meta-modelling process was applied to this data K r but could also have been applied to other quantities of interest for the design, such as the stress-to-deflection ratio or the radial stiffness. This regression problem involves more parameters than the previous application example and allows us to illustrate the interest of the multiplication coefficient form analysis (step 2.1 of the process given in Figure 5). 18), the power coefficient a was constant and equal to 3 for any π i parameters configurations. Thus, only one physical phenomenon drove the evolution of the rotational stiffness: the elasticity of the material. Figure 19, the result of step 2.1, shows that the influence of π 2 and π 3 parameters can be modeled by power laws expressions. The form of the π 1 function cannot be a power law and was modeled here by a polynomial response surface here. The form of the scaling-law-based metamodel was be thus:

ܭ = (ܿ + ܿ ଵ ߨ ଵ + ܿ ଶ ߨ ଵ ଶ )ߨ ଶ ଷ ߨ ଷ ସ ݀ ଷ (29)
The regressions results were (with K r in N.m/rad and d in m) : ܭ = (1.93.10 ଵ + 3.10.10 ଵ ߨ ଵ + 2.30.10 ଵଵ ߨ ଵ ଶ )ߨ ଶ ଶ.ଽ ߨ ଷ ଵ.ହ ݀ ଷ (30) for the first one, where regression was carried out with a and k data (see step 3 section 3.2.3). ܭ = (1.93.10 ଵ + 3.10.10 ଶ ߨ ଵ + 2.30.10 ଵଵ ߨ ଵ ଶ )ߨ ଶ ଷ ߨ ଷ ଵ ݀ ଷ (31) for the second one, carried out on the y data with a least squares fit (see step 3 section 3.2.3).

For the last regression, powers of π 2 and π 3 parameters were set to 3 and 1 respectively. Equation [START_REF] Gogu | Dimensionality Reduction Approach for Response Surface Approximations: Application to Thermal Design[END_REF] presents a simple form that can be handled and implemented in different calculation tools. As shown in Figure 20, the relative error is controlled for a wide range of results with a maximum error of less than 15%. These results are compared to a 3 rd order polynomial function here, which was the metamodel that performed the least poorly of the three alternatives (polynomial response surface, radial basis function, kriging) in the previous application example. The expression of the polynomial response surface required the determination of a large number of parameters: 26 parameters compared to 5 for the proposed form, and showed a maximum error of 21%. The scaling-law-based metamodel allows light and compact models to be obtained, which can be easily integrated into other design algorithms.

CONCLUSION

The design of mechatronic systems requires estimation models to make the link between the design details of each component and the overall design of the system. The article has shown that scaling laws possess interesting properties for performing this role during the sizing tasks. A demonstration of the criteria for obtaining these scaling laws with the Buckingham theorem highlighted the assumptions required for this work. Two of these assumptions, geometric similarity and the uniqueness of the design criteria, limit the generality of their use. A method was proposed to overcome these limitations and generalize the form of scaling laws. The process of regression adapted to this form was described and validated on various examples. The regressions obtained were compared to classical metamodels and showed the interest of the physical basis of scaling-law-based metamodels for preliminary sizing of mechatronic components: their shape is easy to handle while remaining valid over a wide range of sizes, even for prediction or extrapolation purposes. The proposed form of regression can be applied to various fields encountered in mechatronic systems and is illustrated here by simulations on thermal, electromagnetic and mechanical examples.

Figure 1 -

 1 Figure 1 -Mechatronic system (based on reference [3])

1 )

 1 With y the secondary characteristic to be estimated, L the main dimension of the component, and k and a constants. For simplicity of notation, in this article, k means a constant coefficient which may have different numerical values in the different equations.

Figure 3 -

 3 Figure 3 -Brushless motor masses according to the nominal torque.

Figure 4 -

 4 Figure 4 -Bearing masses according to the load capacity.

Figure 5 -

 5 Figure 5 -Meta-modeling process

Figure 6 -

 6 Figure 6 -Power law coefficients calculation

  the possibility of factorizing ݇(ߨ ଵ , ߨ ଶ , … ) by a function of π i . To assess this possibility we introduce the indicator defined by (max(E)-min(E))/(max(log(k)-min(log(k)), where E quantifies the difference between two curves for different values of π i . If the difference between two curves is constant, this indicator is equal to zero and points out the possibility of factorization.• the opportunity of expressing this function by a power law. To assess this, we introduce the indicator (max(s)-min(s))/|(mean(s)|, where s quantifies the slope of a curve. If the slope is constant, the indicator is equal to zero and indicates a power law.

d) Indicator map Figure 7 -

 7 Multiplication coefficient form analysis

Figure 8 -

 8 Figure 8 -Power law coefficients calculation for the test function

Figure 9 -

 9 Figure 9 -Form analysis of multiplication coefficient k for the test function

Figure 10 -

 10 Figure 10 -Sensitivity analysis for the power coefficient a of the test function

a) Effect of noise on step 2. 1 b) Effect of noise on step 3 Figure 11 -

 1311 Figure 11 -Noise effect on indicator map and final regression quality

Figure 13 -

 13 Figure 13 -Mechatronics component from different domains

Figure 15 -

 15 Figure15illustrates the step 2.1 of the process for the cases of winding with or without potting. In the case of winding with potting, the power coefficient varies between 3.42 and 3.49. Without potting, this coefficient varies between 3.23 and 3.48. The 2 cases have a power coefficient varying between 3.5 (demonstrated in section 3.2), characteristic of purely convective thermal transfer, and 3, characteristic of a purely conductive thermal transfer. A coil without potting is naturally more influenced by the conductive effect. For both cases, it was interesting not to consider parameter ܽ(ߨ ଵ ) as a constant.

Range 1 :

 1 L ∈ [20;60] mm and π 1 ∈ [0.5;2.5] ܶ = ݇(ߨ ଵ ܮ) (గ భ ) = 1.35.10 ି ߨ ଵ ଶ.ହଽ ܮ ଷ.ସି.ଵగ భ ା.ଶగ భ మ (T in N.m, L in mm) b) Range 2 : L ∈ [10;100] mm and π 1 ∈ [0.1;5]

2 Figure 16 -

 216 Figure 16 -Quality check of the regressions for different range of study

Figure 17 -

 17 Figure 17 -Geometry of the flexible hinge (2 blades)

Figure 18 -Figure 19 -Figure 20 -

 181920 Figure 18 -Scaling laws with geometrical similarity

Scaling-law-based metamodels for the sizing of mechatronic systems


Frédéric Gueydan
email: frederic.gueydan@um2.fr

Jacques Précigout

Laurent G J Montesi

Strain weakening enables continental plate tectonics

Much debate exists concerning the strength distribution of the continental lithosphere, how it controls lithosphere-scale strain localization and hence enables plate tectonics. No rheological model proposed to date is comprehensive enough to describe both the weakness of plate boundary and rigid-like behaviour of plate interiors. Here we show that the duality of strength of the lithosphere corresponds to different stages of microstructural evolution. Geological constraints on lithospheric strength and large strain numerical experiments reveal that the development of layers containing weak minerals and the onset of grain boundary sliding upon grain size reduction in olivine cause strain localisation and reduce strength in the crust and subcontinental mantle, respectively. The positive feedback between weakening and strain localization leads to the progressive development of weak plate boundaries while plate interiors remain strong.

1/ Introduction

The extrapolation of laboratory flow laws to geological scale suggests a complex layering of brittle and ductile layers within the continental lithosphere [START_REF] Brace | Limits on lithospheric stress imposed by laboratory experiments[END_REF][START_REF] Sawyer | Brittle failure in the Upper Mantle during extension of continental lithosphere[END_REF]. For classical continental geotherm, the upper lithospheric mantle is expected to be brittle and support high stresses. Many analogue and numerical experiments indicate that such a rheological layering is important to reproduce first-order patterns of lithosphere deformation [START_REF] Brun | Narrow rifts versus wide rifts: inferences for mechanics of rifting from laboratory experiments[END_REF][START_REF] Burov | Continental plate collision, P-T-t-z conditions and unstable vs. stable plate dynamics: Insights from thermo-mechanical modelling[END_REF]. In particular, the presence of a brittle uppermost mantle is needed to explain strain localisation at lithospheric scale [START_REF] Buck | Modes of continental lithosphere extension[END_REF][START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF].

However, recent geophysical studies question this classical view of the continental strength layering. Based on earthquakes distribution and elastic thicknesses of the continental lithosphere, including cratons, it has been proposed that the uppermost mantle could behave as ductile instead of brittle [START_REF] Déverchère | Depth distribution of earthquakes in the Baikal rift system and its implications for the rheology of the lithosphere[END_REF][START_REF] Jackson | Strength of the continental lithosphere: time to abandon the jelly sandwich?[END_REF][START_REF] Maggi | A reassessment of focal depth distributions in southern Iran, the Tien Shan and northern India: Do earthquakes really occur in the continental mantle?[END_REF]. However, the mechanical stability of cratons requires that the uppermost mantle supports high stresses [START_REF] Burov | The equivalent elastic thickness (Te), seismicity and the longterm rheology of continental lithosphere: Time to burn-out "crème brûlée"?: Insights from large-scale geodynamic modeling[END_REF]. In addition, the post-seismic displacement field, i.e., the pattern of deformation at the surface of the Earth within weeks to years following an earthquake, suggests that the deep crust is stronger than the lithospheric mantle. [START_REF] Bürgmann | Rheology of the Lower Crust and Upper Mantle: Evidence from Rock Mechanics, Geodesy, and Field Observations[END_REF][START_REF] Thatcher | Temporal evolution of continental lithospheric strnegth in actively deforming regions[END_REF]. Note however that post seismic displacement field may also result from a complex combination of poro-elasticity and fault creep in the seismogenic layer and viscous flow in the lower crust [START_REF] Barbot | A unified continuum representation of post-seismic relaxation mechanisms: semi-analytic models of afterslip, poroelastic rebound and viscoelastic flow[END_REF]. In this case, it may not be used to constrain strength ratios between the ductile crust and lithospheric mantle. Finally, the recent discovery of non-volcanic tremors, i.e., long duration seismic events with small amplitudes, below the San Andreas fault suggests a zone of localized and easily modulated faulting in the lower crust [START_REF] Nadeau | Nonvolcanic Tremor Evolution and the San Simeon and Parkfield, California, Earthquakes[END_REF][START_REF] Thomas | Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas fault[END_REF]. Assuming high temperatures in plate boundaries such as the San Andreas Fault System does not solve this conundrum, as that would imply both a weak mantle and a weak crust.

The above geophysical data question the classical view of continental strength layering and suggest a weak continental lithosphere at plate boundaries. To date, no rheological self-consistent model accounts for both the weakness of plate boundary, which is a prerequisite of plate tectonics, and rigid-like behaviour of plate interiors. A new definition of the continental lithosphere strength is needed to reconcile the apparent contradiction between 1) the mechanical prerequisite of a strong brittle mantle to trigger lithosphere-scale strain localisation, and 2) the low lithosphere strength inferred in actively deforming region, especially in the mantle. Although it may be argued that the presence of fluid and/or shear heating can trigger weak plate boundaries [START_REF] Jackson | Strength of the continental lithosphere: time to abandon the jelly sandwich?[END_REF][START_REF] Thielmann | Shear heating induced lithospheric-scale localization: Does it result in subduction?[END_REF], the existence of long-lived lithosphere-scale inherited weak zone [START_REF] Tommasi | Structural reactivation in plate tectonics controlled by olivine crystal anisotropy[END_REF] suggests that a structural origin for the weakness of plate boundaries is also necessary. Crucially, lithospheric strength must decrease with increasing strain. We present here a quantitative model of strain-dependent lithospheric strength derived from large-strain numerical experiments and guided by field observations that constrain the structural evolution of rocks at various depths in the lithosphere.

2/ Geological constraints

Many processes have been suggested as the origin of strain weakening that triggers strain localization in both the crust and mantle [START_REF] Poirier | Shear localization and shear instability in materials in the ductile field[END_REF][START_REF] White | On mylonites in ductile shear zones[END_REF][START_REF] Montési | A unified description of localization for application to large-scale tectonics[END_REF][START_REF] Regenauer-Lieb | Positive feedback of interacting ductile faults from coupling of equations of state, rheology and thermal mechanics[END_REF]. To date, shear heating has probably received the most extensive theoretical treatment (Brun and Cobbold, 1980;Fleitout and Froidevaux, 1980;Regenauer-Lieb and Yuen, 1998;[START_REF] Kaus | Initiation of localized shear zones in viscoplastic rocks[END_REF][START_REF] Crameri | Parameters that control lithospheric-scale thermal localization on terrestrial planets[END_REF]. However, field observations indicate that metamorphic reactions, fluid infiltration and grain size reduction play a crucial role in enabling shear localization [START_REF] Poirier | Shear localization and shear instability in materials in the ductile field[END_REF][START_REF] Park | Weakening of deforming granitic rocks with layer development at middle crust[END_REF][START_REF] Wintsch | Fluid-rock reaction weakening of fault zones[END_REF]Warren et al., 1996;[START_REF] Gueydan | Analysis of continental midcrustal strain localization induced by reaction-softening and microfracturing[END_REF][START_REF] Drury | Dynamic recrystallization and Strain Softening of Olivine Aggregates in the Laboratory and the Lithosphere[END_REF][START_REF] Gerbi | Magnitude of weakening during crustal-scale shear zone development[END_REF][START_REF] Toy | The role of pyroxenites in formation of shear instabilities in the mantle: evidence from an ultramafic ultramylonite, Twin Sisters Massif, Washington[END_REF][START_REF] Platt | Grainsize evolution in ductile shear zones: Implications for strain localization and the strength of the lithosphere[END_REF]Rennie et al., 2013;[START_REF] Sullivan | Strain localization in homogenous granite near the brittle-ductile transtiion: A case study of the Kellyland fault zone, Maine, USA[END_REF]. It is thus important to characterize the rheological effect of these processes. We will show that geologically-constrained changes in rock texture are able to explain the apparent weakness of the continental lithosphere and the presence of a strong lower crust above a weak mantle at active plate boundaries.

Field and laboratory constraints on weakening processes are summarized here.

In the brittle crust, lubrification of fault zone is achieved by nucleation of new mineral (mica, talc;[START_REF] Holdsworth | PLANETARY SCIENCE: Enhanced: Weak Faults--Rotten Cores[END_REF]. The development of foliated cataclasis enriched in mica leads to strain weakening of fault zone through a decrease of the friction coefficient from 0.6 to 0.1 with strain [START_REF] Collettini | Fault zone fabric and fault weakness[END_REF][START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]. In the ductile crust, the progressive development of layering (shear zone and/or foliation) enriched in mica is also a characteristic of the midcrust and is related to intense weakening the rocks [START_REF] Gueydan | Mechanics of low-angle extensional shearzones at the brittle-ductile transition[END_REF][START_REF] Gueydan | Analysis of continental midcrustal strain localization induced by reaction-softening and microfracturing[END_REF][START_REF] Montési | Fabric development as the key for forming ductile shear zones and enabling plate tectonics[END_REF][START_REF] Wintsch | Fluid-rock reaction weakening of fault zones[END_REF]. Here, we used a two-level mixing theory [START_REF] Montési | A constitutive model for layer development in shear zones near the brittle-ductile transition[END_REF], see Appendix for more details) to calculate the strength of the midcrust as a function of the layering degree, f, which increases from 0 (granitic protolith) to 1 (pure phyllonite) upon deformation. The weakening associated with the development of layering is greatest a low temperature (300°C) while hardening is expected at temperatures higher than 500°C (Fig 1A). This reflects differences in dislocation creep activation energy between mica and quartz/feldspar (Table 1); the mica is much weaker than quartz at low temperature but retains its strength while quartz and feldspar become much weaker than mica at high temperature [START_REF] Kronenberg | Basal slip and mechanical anisotropy of biotite[END_REF]. Furthermore, mica breaks down or melts at high temperature, making the layer development process irrelevant in the lower crust.

Therefore, layering-induced weakening can only occur at midcrustal depth where temperature is below 500°C. Other weakening processes are expected in the deep crust, such as shear heating [START_REF] Thielmann | Shear heating induced lithospheric-scale localization: Does it result in subduction?[END_REF], but are disregard here in order to focus on strain weakening related to microstructural evolution.

In naturally deformed mantle rocks, strain localisation is typically associated with grain size reduction and a switch from grain size-insensitive creep to grain sizesensitive creep [START_REF] Drury | Dynamic recrystallization and Strain Softening of Olivine Aggregates in the Laboratory and the Lithosphere[END_REF]. Recent observations on natural samples highlighted the importance of the grain-size-dependent dislocation-accommodated grain boundary sliding (disGBS) of olivine as controlling the rheology of mantle shear zone [START_REF] Précigout | Strain localisation in the subcontinental mantlea ductile alternative to the brittle mantle[END_REF][START_REF] Warren | Grain size sensitive deformation mechanisms in naturally deformed peridotites[END_REF]. Deformation experiments and numerical investigations have also shown that disGBS could promote strain localisation during dynamic recrystallization [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF][START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]. Indeed, under conditions where disGBS constitutes the dominant deformation mechanism of peridotite, i.e., at temperature lower than 800°C, grain size reduction is associated with a significant drop of strength [START_REF] Précigout | Strain localisation in the subcontinental mantlea ductile alternative to the brittle mantle[END_REF], see Appendix for more details). This amount of strain weakening increases with decreasing temperature and does not occur for temperature larger than 800°C.

(Fig 1B, [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF]).

3/ Large strain numerical experiments

The rheologies and weakening processes described above are implemented in large strain numerical experiments to quantify the relationship between strain and strength at various depths in the lithosphere. Because of the non-linear behaviour of the ductile crust and mantle, numerical modelling is required to capture the interplay between weakening and strain localization and to quantify the strength evolution of the ductile layers of the continental lithosphere. By contrast, strain weakening in the brittle crust simply consists of a change of the friction coefficient from 0.6 to 0.1 due to fabric development [START_REF] Collettini | Fault zone fabric and fault weakness[END_REF][START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF].

3.1/ Numerical results

The numerical experiments follow the evolution of a one-dimensional section of a ductile rock undergoing horizontal simple shear to large strain at a given temperature (Fig. 2). Numerical methods used here for the large strain experiments are very similar to those used in previous studies where techniques for large strain approximation and numerical details of the finite element code SARPP can be founded [START_REF] Gueydan | Mechanics of low-angle extensional shearzones at the brittle-ductile transition[END_REF][START_REF] Gueydan | Continental rifting as a function of lithosphere mantle strength[END_REF]. A week seed of 100 m and with an initial viscosity of 0.9999 times the viscosity of the host rock is imposed in the center of the model in order to initiate strain localisation. Figure 2 shows model resultsextruded to two dimensions to ease readingfor a 100-km wide section of crustal or mantle material sheared by a boundary displacement of 230 km for an overall strain 2.3. In crustal rocks (Fig. 2A), strain localisation occurs at 300°C and 400°C, resulting in a lithospheric-scale shear zone where the local strain reaches 100. The amount of weakening (Fig. 1A) and hence the accumulated strain within the shear zone is maximum at 300°C. Strain localisation results from the progressive development of layering that leads to a decrease in material strength and an increase in strain rate within the shear zone [START_REF] Gueydan | Mechanics of low-angle extensional shearzones at the brittle-ductile transition[END_REF][START_REF] Gueydan | Analysis of continental midcrustal strain localization induced by reaction-softening and microfracturing[END_REF][START_REF] Montési | A constitutive model for layer development in shear zones near the brittle-ductile transition[END_REF]. The crustal shear zones are marked by a degree of layering that progressively tends toward 0.9 with increasing strain while the country rocks strength does not evolve. The shear zone thickness increases with temperature to accommodate the same boundary displacement. For temperature larger than 500°C, the increase in layering leads to strengthening of the crustal rocks and no strain localisation.

A similar behaviour is predicted for the mantle. In this case, viscosity decreases by more than two orders of magnitude at 600°C, leading to the formation of a narrow mantle shear zone (Fig. 2B). Strain localisation results from a weakening induced by the dominance of disGBS during dynamic grain size reduction (Fig 1B and 2B, (Précigout and Gueydan, 2009)). Strain localisation does not develop for mantle temperature larger than 800°C, when disGBS is no longer the dominant deformation process (Fig 1B). In this case, grain size evolution does not induce weakening, and the shear zone viscosity remains close to the initial value.

Figure 3 presents the evolution of the effective viscosity in the developing shear zone in the mid-crust or in the subcontinental mantle. In crustal rocks (Fig. 3A), strain localises at temperatures between 300°C and 400°C. It results from the progressive development of layering that decreases material strength by more than a factor of 10 at strains larger than 1. The asymptotic value corresponds to the complete layering, and hence no further weakening processes can lead to a later drop of strength. At higher temperature, the drop of viscosity is less important and even not present for temperature larger than 500°C. In mantle rocks (Fig 3B), strain localisation induced by the dominance of disGBS leads to a drop of viscosity by more than one order of magnitude at 600°C. The asymptotic value is also reached after strain larger than 1. At this stage, grain size reduction stops and no weakening is expected in the mantle shear zone at larger strain. Strain localisation does not develop for mantle temperature larger than 800°C, when disGBS is no longer the dominant deformation process (Fig 1B). In this case, grain size evolution does not induce weakening, and the shear zone viscosity remains close to the initial value.

3.2/ Strain weakening

Strain localization is associated with a decrease of effective viscosity from an initial value  0 to an asymptotic value  ∞ in ductile shear zones or a reduction of strength from  0 to  ∞ in the brittle crust. The effective viscosity is defined by the ratio between the shear stress and the overall strain rate that is here 10 -15 s -1 . Figure 4 summarizes these results by presenting a measure of strain weakening, , as a function of depth, where =1-( ∞ / 0 ) or =1-( ∞ / 0 ) in brittle or ductile layers of the crust. For the upper crust, a change from 0.6 to 0.1 implies a strain weakening of 83% related to the progressive development of foliated cataclasites (Fig. 4).

Weakening is even more intense in the middle crust where, for example, viscosity decreases from 2x10 23 Pa.s to 10 22 Pa.s and implies a strain weakening of 95% (Fig 4). These relative high values of viscosity compared to what is inferred in most plate boundaries [START_REF] Thatcher | Temporal evolution of continental lithospheric strnegth in actively deforming regions[END_REF]Pollitz, 2008, Bürgmann and[START_REF] Bürgmann | Rheology of the Lower Crust and Upper Mantle: Evidence from Rock Mechanics, Geodesy, and Field Observations[END_REF] are due to the selected value of mean strain rate that would however not change the quantification of strain weakening, since it is a ratio of viscosities. Strain weakening remains close to 90% in the crust at temperature lower than 400°C and then progressively decreases to vanish at 500°C. No strain weakening is expected in the deep crust, since shear heating is disregarded. In the lithospheric mantle, strain weakening reaches again 90% for temperature lower than 700°C and then decreases to vanish at 800°C [START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF].

The evolution of strength predicted by our numerical models as a function of strain can be approximated by an exponential relationship as follow:

⁄ (eq. 1)
in which is the strain and is a parameter that represents the characteristic strain over which the fabric of deforming rocks changes according layering development or grain size reduction. The large strain experiments (Figs. 23) are represented well with this general form and with a critical strain of ~ 0.5 for both the ductile crust and the ductile mantle. For seek of simplicity, the critical strain for the brittle crust is also set to 0.5. However, note that the critical strain depends on the poorly constrained kinetics of metamorphic reaction for the crust and of grain size reduction for the mantle that remain to be quantified by large strain experimental deformation in laboratory.

Combining eq 1 with the definition of the strain weakening factor  yields an unified description of the lithosphere strength that becomes a function of strain , temperature T, and strain rate ˙ e :

̇ ̇ ⁄ ̇ (eq. 2)
The initial strength ̇ is defined in the ductile layers for the protolith (granitic composition or an olivine aggregates in the crust or the mantle, respectively) according to the classical strain-rate and temperature sensitive flow law such as dislocation creep (Appendix section). In the brittle upper crust, the initial strength obeys the Mohr-coulomb failure criterion. It is insensitive to temperature and strain rate but depends on confining pressure, and therefore depth. The values of the strain weakening coefficient  are those shown on Figure 4 and therefore are a function of depth.

4/ Continental lithosphere strength profiles during strain localization

The unified model of strength as a function of strain, strain rate and temperature (Eq. 2) is now used to characterize the evolution of the continental lithosphere strength during progressive strain localization. The development of plate boundary is represented by an increase in strain rate from 10 -17 s -1 to 10 -13 s -1 coeval with a strain increase from 0 to 2 (Fig. 5). The continental geotherm is such that the Moho temperature is 600°C. We ignore the possibility of shear heating to focus on microstructural changes that accompany localisation.

At low strain rate and strain (Fig. 5A), the continental strength profile shows a classical mechanical layering, with two high strength layers: the upper crust and the uppermost mantle. No brittle mantle is predicted at such a low strain rate. As strain rate increases to 10 -15 s -1 and strain reaches 1, two contrasting features appear (Fig. 5B). Strength decreases in the upper crust, mid-crust and the uppermost mantle where strain weakening is expected because of brittle faulting, layering development and grain size reduction coupled to disGBS, respectively (Fig. 5B). By contrast strength increases markedly in the deep crust where no strain weakening is suspected. The same features are strongly enhanced at the final strain rate and strain (10 -13 s -1 , ε=2, Fig. 5C). Overall, strain weakening is associated with a decrease in lithosphere strength as strain rate increases to counteract the strength increase otherwise produced by the fundamental strain-rate hardening characteristics of ductile flow laws.

The positive feedback between strain localization and weakening leads to a reversal in strength stratification of the lithosphere. The lower crust, which is initially the weakest layer, becomes the strongest, while the upper mantle and upper to middle crust become the weakest levels. Thus it is possible to reconcile the apparently contradictory evidence for a high strength upper mantle needed to develop plate boundaries and the evidence of a strong lower crust in active deformation zones.

5/ Discussion and conclusion

The strength model developed above shows that actively deforming regions such as plate boundaries, are characterized by a weak upper-mid crust and weak lithosphere mantle surrounding a high strength deep crust (Fig. 5C), where the zones of weakness are related to microstructural evolution. The weakness of continental plate boundaries is consistent with geophysical observations that suggest very low elastic thickness in these regions [START_REF] Audet | Dominant role of tectonic inheritance in supercontinent cycles[END_REF].

Strain weakening in the ductile middle crust permits mechanical decoupling between lower crustal flow and upper crust faulting, which is essential to promote localised faulting and large displacement along faults [START_REF] Lavier | Half graben versus large-offset low-angle normal fault: importance of keeping cool during normal faulting[END_REF][START_REF] Huismans | Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening[END_REF][START_REF] Nagel | Symmetric alternative to asymmetric rifting models[END_REF][START_REF] Schueller | Brittle-ductile coupling: Role of ductile viscosity on brittle fracturing[END_REF][START_REF] Schueller | Mechanics of the transition from localized to distributed fracturing in layered brittle-ductile systems[END_REF]. The absence of weakening in the deep crust yields to a progressive increase of its strength during strain localization so that it becomes the strongest level of the lithosphere, as is indicated by post-seismic creep measurements in California [START_REF] Bürgmann | Rheology of the Lower Crust and Upper Mantle: Evidence from Rock Mechanics, Geodesy, and Field Observations[END_REF][START_REF] Thatcher | Temporal evolution of continental lithospheric strnegth in actively deforming regions[END_REF].

The long-term strength of the continental lithosphere described here may be further modified by interaction with the seismic cycle. Figure 6 proposes a structural interpretation of the character of seismicity at a continental strike slip fault like the San Andreas Fault. Major earthquakes take place in the seismogenic upper crust.

Immediately below that level, the localizing midcrust, which acts as a decoupling layer, constitutes a region where fluid-rock interaction is dominant, allowing metamorphic reaction, development of layering, and related weakening [START_REF] Holdsworth | PLANETARY SCIENCE: Enhanced: Weak Faults--Rotten Cores[END_REF]. In this layer, weakening would generate relatively narrow ductile shear zones where transient changes in pore fluid pressure and shear stress may induce microseismicity [START_REF] Rolandone | The evolution of the seismicaseismic transition during the earthquake cycle: Constraints from the timedependent depth distribution of aftershocks[END_REF][START_REF] Gratier | Deformation in the presence of fluids and mineral reactions: Effect of fracturing and fluid-rocks interaction on seismic cycles[END_REF].

Post-seismic deformation following large earthquakes is detected at large distance from the main shock, implying the presence of a strong layer at depth most likely in the lower crust [START_REF] Freed | Far-reaching transient motions after Mojave earthquakes require broad mantle flow beneath a strong crust[END_REF][START_REF] Lindsey | Geodetic slip rates in the southern San Andreas Fault system: Effects of elastic heterogeneity and fault geometry[END_REF]. Non-volcanic tremors at 20 to 40 km depth below the San Andreas fault [START_REF] Nadeau | Nonvolcanic Tremors Deep Beneath the San Andreas Fault[END_REF][START_REF] Thomas | Tidal triggering of low frequency earthquakes near Parkfield, California: Implications for fault mechanics within the brittle-ductile transition[END_REF] suggests that shear failure is possible within the dominantly creeping deep crust. The easy modulation of tremors activity in response to different loading (Shelly andHardebeck, 2010, Thomas et al., 2012) furthermore implies that the lower crust is very weak and close to failure. Transient loading during earthquakes and slow slip events in the mid-crust and in the upper mantle may increase pore fluid pressure and lead to shear failure in the otherwise ductile lower crust, generating tremor activity (Fig. 6). Further studies, combining geological, mechanical, and seismological approaches, are needed to unravel the precise mechanics of transient loading, fluid-rock interaction, pore fluid pressure and shear failure.

Finally, the mantle opposes little resistance to deformation as evidenced by the absence of earthquakes in the mantle and also postseismic displacements pattern. Our rheological model supports the idea of a broad mantle shear zone, probably 50 km wide below the fault zone [START_REF] Freed | Far-reaching transient motions after Mojave earthquakes require broad mantle flow beneath a strong crust[END_REF]Fig 6). The exact thickness of this mantle shear zone would depend on the degree of heterogeneity [START_REF] Vauchez | Faults (shear zones) in the Earth's mantle[END_REF], the amount of coupling between crustal deformation and mantle deformation [START_REF] Schueller | Mechanics of the transition from localized to distributed fracturing in layered brittle-ductile systems[END_REF] and the amount of weakening within the lithospheric mantle. The rheological model developed here provides the foundation for new study of shear zone development that would unravel the expected structure of continental fault zones at depth.

6/ Acknowledgements

We thank Roland Bürgmann, an anonymous reviewer, and editors Evgenii Burov and Laurent Jolivet for their constructive comments on this manuscript.

LGJM was supported by grant NASA PGG NNX10AG41G.

7/ Appendix: Strain-dependent rheologies

6.1/ Midcrust rheology

The crustal protolith is assumed to have a granite-like mineral assemblage, constituted of 40% quartz, 50% feldspar, and 10% mica homogenously distributed throughout the rock. The strength of the protolith,  protolith , is the weighted sum of the strength of its constitutive minerals assuming that geometrical incompatibilities in a homogeneous rock force all the phases to adopt the same strain rate:

 protolith = 0.4 quartz + 0.5  feldspar + 0.1  mica (eq. 3)

where  quartz ,  feldspar and  mica are defined by the experimentally calibrated flow such as the dislocation creep ̇ (eq. 4)

where A, Q, and n are the pre-exponential constant, the activation energy and the stress exponent and depend on the rock type (Table 1).

After large deformation, the protolith becomes a mylonite, with higher mica content due to fluid-rock interactions and a strain-induced layered structure. If reaction reaches completion, the strength of the shear zone would be given by the rheology of mica:

 shear zone =  mica (eq. 5)

At intermediate stages, a two-level mixing theory [START_REF] Montési | A constitutive model for layer development in shear zones near the brittle-ductile transition[END_REF] is used to estimate the strength of the aggregate as a function of the degrees of layering f:

=(1-f) protolith +f shear zone (eq. 6)

Although a pure mylonite would have f=1, we consider the final stage of evolution to be f=0.9 to avoid strong numerical instabilities that arise when weakening is too intense.

/ Mantle rheology

Based on field observations in the Ronda peridotite (southern Spain) and the experimental data of [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF], [START_REF] Précigout | Strain localisation in the subcontinental mantlea ductile alternative to the brittle mantle[END_REF] proposed a strain-dependent mantle rheology in which the overall strain rate ̇ is the sum of four ductile deformation mechanisms: dislocation creep ( ̇), diffusion creep ( ̇), disGBS ( ̇) and exponential creep ( ̇) [START_REF] Hirth | Rheology of the Upper Mantle and the Mantle Wedge: A view from the Experimentalists[END_REF][START_REF] Drury | Dynamic recrystallization and Strain Softening of Olivine Aggregates in the Laboratory and the Lithosphere[END_REF][START_REF] Précigout | Strain localisation in the subcontinental mantlea ductile alternative to the brittle mantle[END_REF]. Each mechanism contributes to the bulk strain rate ( ) of an olivine aggregate according to:

̇ ̇ ̇ ̇ ̇ (eq. 7)
where the corresponding flow laws are: ̇ (eq. 8)

̇ (eq. 9)

̇ (eq. 10)

and ̇ (eq. 11)

In the above equation,  is the stress [MPa], d is the grain size

[µm], T is temperature [K], R

References e

Figure 1 :

 1 Figure 1: Field observation in naturally deformed rocks allowing the definition of

Figure 2 :

 2 Figure 2: 1D large strain simple shear numerical experiments, shown in 2D to ease of

Figure 3 :

 3 Figure 3: Evolution of effective viscosity during large strain experiment (inset and Fig.

Figure 4 :

 4 Figure 4: Strain weakening quantified by large strain experiments (Figs 2-3) as a

Figure 5 :

 5 Figure 5: Three strength profiles for low (A), medium (B) and high (C) strain rate and

Figure 6 :

 6 Figure 6: Schematic diagram of the mechanical layer of a continental strike-slip plate

  Grainsize reduction can only occur within grain size/stress domains where dislocation creep is effective enough, i.e., dislocation creep and disGBS, up to the boundary of the diffusion creep field, where we assume that only grain growth can occur.

is the gas constant, and A, Q and n are the material parameters indicating the pre-exponential constant, the activation energy and the stress exponent for each mechanism identified by their respective indexes (r=Dislocation creep; d=Diffusion creep, g=DisGBS creep and e=exponential creep). Finally, m d , m g , and  are respectively the grain size exponent for diffusion creep and disGBS, and a constant parameter for the exponential creep (Peierls stress;

[START_REF] Goetze | The mechanisms of creep in olivine[END_REF]

(Table

1

).

 [START_REF] Wilks | Rheology of some continental lower crustal rocks[END_REF]


J.-F Colomer

C Stephan

S Lefrant

G Van Tendeloo

I Willems

Z Konya

A Fonseca

Ch Laurent

J B Nagy

Large-scale synthesis of single-wall carbon nanotubes by catalytic ž / chemical vapor deposition CCVD method

´Ábstract Ž

Introduction

In addition to the laser evaporation and electric arc discharge techniques, the catalytic process is a very efficient method to produce multi-wall carbon w x nanotubes 1,2 . The first two techniques allow the w x production of SWNTs in high quantities 3,4 . Many w x previous reports 5-10 have shown that the catalytic way could be a possibility to produce nanotubes on a large scale at a very low cost.

w x Biro et al. 5 reported the presence of SWNTs by Ž . scanning tunneling microscopy STM during the ) Corresponding author. Fax: q32-81-72-46-00; e-mail: jean-francois@fundp.ac.be synthesis of multi-wall nanotubes by catalytic decomposition of hydrocarbons.

w x Peigney et al. 6 also reported the synthesis of a mixture of single-and multi-wall nanotubes by decomposition of H rCH on Fe alumina nanocom- . mina with varying efficiency, depending on the system used. The abundance of SWNTs seems to increase when binary alloy nanoparticles, for the same metals, are used as in the case described by w x Flahaut et al. 11 . The main inconvenient of this method is to estimate the yield of the reaction in SWNTs because of the presence of catalyst and the Ž other forms of carbon multi-wall nanotubes Ž . MWNTs , amorphous carbon, carbon nano-. spheres, . . . . As a conclusion, in all cases described in the literature, the synthesis of SWNTs seems to be linked to the very small size of the catalyst particles. A mechanism for the nucleation of both single-and multi-wall nanotubes obtained by CCVD from catw x alytic particles has already been reported 7,10 .

Note that pyrolisis of hydrocarbons with metallocenes has also been successfully applied for the w x production of SWNTs 12,13 in large quantities, but this way does not use the metal-supported catalysts, like the others described above. Ž . produce both single-and multi-double-wall carbon nanotubes in a 1:1 ratio. The carbon nanotubes were extracted by dissolution of MgO and part of the Ž Co catalyst in the HCl aqueous solution 36% room . temperature . Indeed, MgO presents the advantage over other supports such as alumina, spinel and silica that it can be readily dissolved in acids. It was w x reported 14,15 that the mild acid treatment does not damage the carbon nanotubes, more than 80% of which are either SWNTs or DWNTs.

In the present work, we describe the large-scale synthesis of SWNTs over supported catalyst MetalrMgO where Metal is Co, Ni, Fe and a binary mixture Co-Fe. The possibility to obtain a sample Ž . containing high yields of SWNTs 70-80% , by this way is an improvement. This fact is explained by the synthesis conditions, which allow the formation of Ž SWNTs in abundance which some double-wall nan-. otubes , and by the use of the support such as MgO known to be easily removed by an acidic treatment Ž . HCl . Using the transmission electron microscopy Ž . TEM , the quality of the SWNTs and the nature of contaminating phase were checked. The SWNTs were also characterized by Raman spectroscopy.

Experimental

The different catalysts are prepared by impregna-Ž . tion of MgO Vel with ethanol solution of metal Ž . Ž . salts Co, Ni, Fe or a mixture of metal salts Co-Fe in the appropriate concentration. The material is sonicated during one hour and the ethanol is removed via a rotary evaporator. Then, the material is dried at 1308C during 12-15 h and ground into a fine powder. Finally, a portion of each catalyst was Ž y1 hydrogenated in H rN H flow: 30 ml min , N 2 2 2 2 y1 . flow: 80 ml min .

The synthesis of SWNTs was carried out in a Ž fixed-bed reactor quartz tube of 50 mm diameter, 80

. cm in length in a carbolite horizontal reactor at 10008C, with a typical reaction time of 10 min. For each synthesis, a quartz boat containing about 1 g of catalyst was placed in the center of the reactor. The carbon source was methane with hydrogen as carrier gas, and the flows of the mixture methanerhydrogen were 75r300 ml min y1 . After the reaction, the Ž . catalyst support and metal particles is removed by acidic treatment in concentrated hydrochloric acid.

Ž . Typically, the crude sample about 1 g is placed in 50 ml of concentrated HCl under sonication for 15 Ž min. It is filtered through an HPLC filter FP-Vericel . membrane filter, porosity 0.2 mm washed with distilled water to fix the pH at 6-7 and finally, dried at 1308C.

After this purification, the carbon yield can be calculated as follows: Carbon yield % s 100 m rm , Ž .

Ž .

pure crude where m is the weight of the sample after procrude duction of SWNTs, and m the weight of the pure sample after acidic treatment. Finally, the sample for TEM was prepared by sonication of about 3 mg of material in ethanol and a Raman experiments have been performed with a Raman spectrophotometer Jobin Yvon T64000 in ambient atmosphere and at room temperature. The Ž spectra are recorded using a red excitation line 676.4 . nm .

Results and discussion

Observations by low-resolution microscopy Ž . TEM: Philips CM 20 show entangled carbon filaments and the abundance is very high in all the Ž . observed samples Fig. 1 . From several TEM images, the yield of these filaments can be estimated approximately to be of the order of 70-80%, depending on the catalyst used. High-resolution trans-Ž mission electron microscopy HREM using a JEOL . 200 CX images of the sample show that each filament consists of bundles of SWNTs but confirm also the presence of isolated SWNTs in a large Ž . quantity Figs. 2 and3 . One of these SWNTs bundles, forming crystallite-like entities organized in a two-dimensional lattice, is shown in Fig. 3b.

Moreover, the diameter of the SWNT can be measured by transmission electron microscopy. Typically, the diameters of isolated SWNTs are 1-5 nm. For the SWNTs aligned in the bundles, the diameter values vary between 0.8 and 2 nm, and the estimated diameters distribution is larger than those synthe- sised by arc discharge or laser evaporation. The SWNTs diameters were also studied in detail using Ž . Raman spectroscopy Fig. 4 . The spectrum exhibits unambiguously the characteristic frequencies of Ž Ž SWNTs, but the corresponding bands low 100-250 y1 . Ž y1 . . cm and high 1500-1600 cm frequencies are reduced by the presence of other carbon forms like Ž multi-wall nanotubes or amorphous carbon 1337 y1 . cm , band D . The high-frequency bands can be decomposed in one main peak around 1581 cm y1 , with a shoulder around 1540 cm y1 . This shoulder is more important in the case of SWNTs produced by arc discharge, where a separate peak can be observed w x 4 . The high-resolution spectrum obtained in the low-frequency domain shows several components at 121, 133, 150, 155, 172, 175, 190, 196, 211, 221 cm y1 . The spectrum in this frequency domain is very sensitive to the tube diameters, where the fre-Ž . quency increases with decreasing tube diameter d . Moreover, these bands are due to A symmetry g y1 Ž .

Ž . Ž . Ž .
w x mode where n cm s 2238rd A 16 and the results reflect a distribution in diameters. Then, SWNTs with small diameters can be observed with the characteristic bands 211, 221 cm y1 , where the diameters 10.6 and 10.1 A are calculated from the formula above.

The so characterized sample, by Raman spectroscopy and by HREM, shows also the presence of Ž small amounts of impurities like MWNTs which are . often double-wall nanotubes , metal particles encapsulated in the nanotubes or in polyhedral graphitic nanoparticles and amorphous carbon and seems to be less pure than those produced by the other techniques.

In the sample, the SWNTs tips are observed for isolated SWNTs and catalyst particles can be present Ž . or not in the tube tip Fig. 5 . This fact is in agreement with the suggested mechanism, where the growth mechanism of isolated SWNTs is close to that of the MWNTs in the CCVD process, and the diameter of the SWNT is determined by the size of the catalytic metal particle attached to their end w x 7,10 . The large diameter distribution of isolated SWNTs is also explained by this fact.

The dispersion of metal on the support surface influences the metal particle size and hence the presence of isolated SWNTs or SWNTs bundles. The Ž . bigger metal particles ) 15-20 nm are not active in the synthesis of SWNTs and are often found encapsulated by several graphitic layers. All the other metal particles, having a size smaller than that of the big particles are active in the synthesis of SWNTs bundles. The hypothesis of SWNTs bundles growth is that each SWNTs bundle grow from a single metal particle.

Concerning the catalysts, the nature of support Ž . commercial magnesium oxide is very important because it is easily removed and this fact allows one to estimate the yield of the reaction and to characterize with more precision the nature of carbon deposit: quality of SWNTs, nature of impurities, etc. It is the main advantage of this support with respect to the other supports used in the synthesis of SWNTS like alumina or silica.

Ž . For the catalysts used Table 1 , the activity seems be dependent on the nature of metal, but in all cases Ž and for the two different preparations hydrogenated . or not , the yield is around 7%, except for the Ni Ž . catalysts where the yield is very low -1% . The best activity found is for the catalyst prepared with Ž . the Co-Fe mixture. But in this case presence of Fe and so for the Fe catalysts, the impurities like MWNTs and encapsulated metal particles seem to be more abundant than with the Co catalysts. Moreover, Ž when the weight percent of metal increases from 2.5 multi-wall carbon nanotubes could be produced when the temperature is between 800 and 10008C.

Conclusion

Ž We have shown that the large-scale production 1 y1 . Ž . g day of single-wall carbon nanotubes SWNTs is possible by CCVD method similar to the two other techniques used until now, arc discharge and laser evaporation. Indeed, the method using the catalytic decomposition of methane over well-dispersed metal particles supported on MgO at 10008C has been adapted to produce large quantities of SWNTs. The synthesis yield was easily estimated, after removing the support by a simple acidic treatment to obtain a product containing SWNTs of high purity. The SWNTs content in the purified samples was estimated approximately to be of the order of 70-80%.

For a 10 min reaction, we obtain around 50 mg of carbon deposit and for 20 reactions day y1 , the daily production is of about 1 g of product containing essentially SWNTs. This method of synthesis could replace the other two techniques of SWNTs produc-Ž . tion, because of the lower temperature used 10008C , of the low cost of production and of its industrial application potential.

In the near future, the purification of SWNTs produced by this way will be optimised to obtain pure SWNTs in large quantities.

  Finally, the production of SWNTs in-Ž . cluding double-wall nanotubes DWNTs on Mo and w x Mo-Fe alloy has also been reported 9 . By this method, the choice of the support and the metals or metal mixtures seems be determinant in the SWNTs w x synthesis. Colomer et al. 10 have shown the Ž presence of SWNTs for different metal pairs Fe, Co, . Ž Ni, binary mixture, . . . rsupport silica and alu-

Fig. 1 .

 1 Fig. 1. Low-resolution image of abundant carbon filaments pro-Ž . duced on Co 2.5% wt rMgO catalyst.

  Fig. 2. A and B . High-resolution images of carbon filaments produced on Co 2.5% wt rMgO catalyst showing that these filaments are ropes of SWNTs. Higher-resolution images of A SWNTs organized in bundles or isolated synthetised on Co 2.5% wt rMgO catalyst, and B cross-section of SWNTs bundle showing the triangular lattice arrangement. Ž Ž . . Fig. 4. Raman spectrum of SWNTs on Co 2.5% wt rMgO catalyst recorded at room temperature, using an excitation wavelength of 676.4 Ž y1 . nm. Inset: the low-frequency range 50-250 cm has been given in detail.

Fig. 5 .

 5 Fig. 5. High-resolution image of single-wall carbon nanotube with a small catalytic particle encapsulated to the tip. Double-wall Ž carbon nanotube can be observed on the left of the picture this Ž . . sample was produced on Co 2.5% wt rMgO catalyst .

Table 1

 1 Yields of SWNTs synthesis depending on the catalysts used

	Metals	Weight	Yields for	Yields for
	used	percent	as prepared	hydrogenated
		of metal	catalyst	catalyst
	Co	2.5%	5.5%	5.7%
		5%	6.5%	4.7%
	Fe	2.5%	6%	5.9%
		5%	5.8%	5%
	Co-Fe	2.5-2.5%	7.6%	8.3%
	. to 5% , the yields are similar, but the quality of the
	samples seems to be different, because the increase
	of the metal content in the catalyst increases the
	encapsulation of metal in the graphitic nanoparticles
	and hence the quantity of SWNTs is smaller.
	The conditions of SWNTs production, i.e.,
	methane used as hydrocarbon, hydrogen as carrier
	Ž gaz, reaction time 10 min and temperature 10008C . Ž .
	allow the formation of SWNTs containing only a
	small amount of amorphous carbon and polyhedral
	graphitic nanoparticles. In fact, the reason to produce
	mainly SWNTs in abundance in the present work is
	Ž the constant reaction temperature 10008C during .
	the reaction. It is possible that under other reaction
	w conditions 14,15 , i.e., increasing until 10008C and x Ž .
	decreasing the reaction temperature, under H rCH , 2 4

Acknowledgements

The authors are grateful to the European Commission, TMR contract NAMITECH, ERBFMRX-Ž . CT96-0067 DG12-MIHT and thank the Belgian Programme on Inter University Poles of Attraction initiated by the Belgian State, Prime Minister's Of-Ž . fice of Science Policy Programming 4r10 .

w x 1 M.J. Yacaman, M.M. Yoshida, L. Rendon, J.G. Santiesteban, Ž . Appl. Phys. Lett. 62 1993