J Labroqu

R Duvigneau
email: regis.duvigneau@inria.fr

T Kloczko

J Wintz

INTERACTIVE COMPUTATION AND VISUALIZATION TOWARDS A VIRTUAL WIND TUNNEL

INTRODUCTION

Computational Fluid Dynamics (CFD) software is now commonly used for aerodynamic studies as a replacement of traditional experimental facilities. CFD codes take benefit from the increase of the computational facilities and the maturity of numerical methods. However, they are usually based on old-fashioned software architectures, inherited from the 80s. They are typically built as static codes, independently from software components that are used to construct the geometry and the grid, or visualize flow fields. As a consequence, the whole simulation process is usually quite complex, including several different phases, and is restricted to expert users. Moreover, it is not possible for the user to interact with the computed flow, which makes a significant difference with an experimental wind tunnel. In the latter case, the user can for instance adjust inflow parameters, or modify some geometrical characteristics such as wing incidence, while interactively observing the flow evolution. The studies on interactive CFD simulation are quite uncommon, and a very few papers can be found on this topic [START_REF] Kreylos | Interactive visualization and steering of cfd simulations[END_REF][START_REF] Schirski | Vista flowlib a framework for interactive visualization and exploration of unsteady flows in virtual environments[END_REF][START_REF] Wenisch | Collaborative and interactive cfd simulation using high performance computers[END_REF].

In this study, we investigate the use of a modern software architecture in the context of CFD, which allows the user to interact with the computation, by modifying physical or numerical parameters during the computation and visualize the impact on the flow. Our objective is to evaluate the interest of such a software architecture and measure the possible benefit for scientific studies.

We describe in a first section the proposed software architecture and its new features, that allow the user to interact with the computation and visualization. Then, the application to compressible flow simulation is considered. In particular, we show that various implementations can be envisaged, with different advantages and drawbacks. The possible interactions with the computation and visualization are described. In a third section, some illustrations are proposed. Finally, the bene-fit of the use of these new features are discussed from the point of view of the CFD practitioner.

SOFTWARE ARCHITECTURE

Platform and plugins system

Traditionally, CFD codes are only composed of a main program built with some static libraries. The main program describes more or less the workflow that corresponds to the simulation, whereas the librairies contain the functions necessary for each task of the workflow. This straightforward program architecture has been set up during the 70s and 80s and has not really be questioned so far, although the complexity of the computations has grown significantly.

However, it appears clearly that this architecture has several drawbacks. First, the achievement of numerical experiments requires modularity, which allows to easily test different numerical methods (e.g. Roe flux v.s. HLLC flux), or modify physical parameters (e.g. inlet velocity profile). This ability to change some parts of the code is critical to precisely benchmark numerical methods or models. Obviously, the traditional static architecture is not well suited to this purpose: usually, the implementation of a new method foo2 yields the addition of new functions and variables, that will juxtapose the existing ones foo0, foo1. Then, the code is always growing and becoming less and less readable.

Developments in CFD are now mainly collaborative projects, because it is more and more difficult for one isolated person to master all modeling, numerical and computational aspects related to complex simulations. Unfortunately, the traditional architecture is a real burden for collaborative development: since these codes are always growing (as explained above), the exchange of pieces of code and the maintenance is time consuming for long term and large scale projects. Typically, for each new method foo2, someone should verify if the proposed implementation is correct and will not generate conflicts with all other existing methods foo0, foo1 (for instance in case of addition of an argument).

Finally, the complexity of these codes make the introduction of a new user tedious and time consuming: to modify a unique method, a new user should usually understand a large part of the code, even if it is not of interest for his study.

To overcome these limitations, we have initiated the NUM3SIS project (http://num3sis.inria. fr), whose software architecture is based on the distinction between the platform and plugins. The platform, written in C++ language, gathers all components dedicated to numerical simulation, that are considered as common to various computations. In practice, it consists of a set of abstractions, that can represent data or processes, commonly used in simulation. For instance, the core of the platform is composed of abstractions for grids, fields, flux computations, finite-elements, etc. As abstraction, the numerical methods related to these objects are not implemented in the platform.

On the contrary, a plugin contains a possible implementation of an abstraction defined in the platform. Note that a library written in a different language can be embedded into a plugin. In practice, plugins are dynamic libraries used by the platform at runtime. For instance, foo0, foo1 and foo2 can be three different plugins (possibly based on existing libraries) implementing the abstraction of the method foo, defined in the platform (see Fig. 1). This approach has several advantages: first, all the methods are not aggregated in a unique code, which improves readability and modularity. Then, it proposes a nicer framework for collaborative development, because development of new plugins can be conducted independently from the platform or other plugins. Moreover, templates can be proposed to speed-up the coding phase of plugins. Finally, an easy benchmarking procedure can be carried out by implementing the methods to be compared into different and independent plugins. One should underline that it is possible to change a parameter of a plugin, or the plugin itself, at runtime. This modifies the forthcoming computation and allows the user to interact with the simulation. This point will be detailed latter, with some examples.

Visual programming

As already mentioned, it is not easy for a new user to implement a method in an existing CFD code. This is particularly the case for people who are more familiar with mathematics than programming. This is dommageable because it is a real obstacle to the improvement of numerical methodologies. Therefore, we have introduced a visual programming tool in the platform, that allows to build a high-level computational scenario without writing any line of code. Thus, for each abstraction defined in the platform, a visual wrapping (node) is introduced, which can be handled in the composer space, and connected to other nodes to construct the desired computational scenario (see Fig. 2). For a given node, several plugins may exist, that correspond to various implementations, and can be selected using the composer. This approach is hierarchical, in the sense that a node can contain itself a composition based on other nodes. Moreover, particular nodes represent control structures, such for loops, if then conditions, etc. This visual programming approach really facilitates the introduction of a new user. Moreover, it represents a valuable tool to prototype a new computational scenario. It should be underlined that the frontier between visual programming (using the composer) and inline programming (implementation in plugins) completely depends on the choice of the user: it is possible to program entirely the simulation process using the composer and, on the contrary, it is also possible to embed a whole simulation process into a unique plugin. Of course, an intermediate choice is usually more useful. Performance studies have been conducted to quantify a possible loss of performance due to the use of such a visual programming tool. However, it has been found that this loss of performance is negligible.

Graphical interface and scientific visualization

The graphical interface is composed of four main spaces: the composer space that allows to visually build a computational scenario, the visualization space used for scientific visualization of computational results and two other spaces dedicated to data exchange and parallel computing.

The composer space permits the user to select a node among a set of existing nodes, drop it in a working space, select the desired implementation in a corresponding plugin (see Fig. 2), and connect it with other nodes to construct a given computational scenario.

The visualization tool is based on the wellknown VTK library (http:://www.vtk.org) and contains all usual features to visualize scalar or vectorial fields, glyphes or streamlines, grid nodes, faces, edges, etc (see Fig. 3). The main important feature is that an object can be visualized as soon as it is computed. Actually, the view itself is considered as a node in the composition which is fed by grids, fields, etc. Therefore, each time a field is updated, the modifications can be seen in the visualization space. This feature, associated with the possibility to change a parameter or a plugin at runtime, makes the interactive computation and visualization possible.

Stereoscopic visualization

The researches related to methods used to visualize 3D fields have known a growing interest for the last years. Indeed, simulated flows are more and more complex and the use of sophisticated turbulence models, such as LES (Large Eddy Simulation), DES (Detached Eddy Simulation) or VMS (Variational Multi-Scale) models requires visualization tools that help to capture the main char-acteristics of the flow. Among all the possible approaches, the use of virtual reality facilities is explored. INRIA Sophia Antipolis M éditerran ée center has an immersive space facility, which regroups two display devices, a CadWall for its ease of use and an iCube for its immersion quality. This facility is dedicated to research in virtual reality, but can also be used by non-specialists to explore the possibilities of these new visualization devices.

The CadWall employed for the present experiments can be seen as a single screen of 3528x1200 pixels. Two images are simultaneously projected onto the wall and allows to create a 3D perception of the objects by using specific glasses. A dedicated software layer has been introduced in the platform in order to make the generation of such a stereoscopic view possible. Then, the platform can be used with the CadWall, or with any classical screen, without any modification.

APPLICATION TO COMPRESSIBLE FLOW SIMULATION

We explain in this section how the proposed software architecture is used for compressible flow simulation. We present first the numerical methods employed and then how they are implemented in the platform and plugins. Note that the platform is not devoted to CFD, fluid mechanics is only an application among others. At the present time, the platform is used by two INRIA Project-Teams for computations in aerodynamics (finite-volume method), electromagnetics (discontinuous Galerkin method), pedestrian traffic modeling (finite-volume method) and thermal conduction (finite-element method). Moreover, the computational scenario is not restricted to simulation, but other computations can be carried out, on the basis of the same software components.

Modelling and numerics

Modeling The equation solved for compressible flows simulation are the Navier-Stokes equations. In conservative form, they can be written as:

∂ t W c + ∇ • F (W c ) = ∇ • N (W c ) + S(W c ) (1)
with the conservative variables W c = (ρ, ρu, ρv, ρw, ρE) T , the inviscid flux

F (W c ) = (F (W c ), G(W c ), H(W c )) T , the viscous flux N (W c ) = (R(W c ), S(W c ), T (W c )) T and the source term S(W c ).
The components of the inviscid flux in the global frame R 0 ( x, y, z) are:

F (W c ) =       ρu ρu 2 + p ρuv ρuw u(ρE + p)       , G(W c ) =       ρv ρvu ρv 2 + p ρvw v(ρE + p)       H(W c ) =       ρw ρwu ρwv ρw 2 + p w(ρE + p)      
The components of the viscous flux are:

R(W c ) =       0 τ xx τ xy τ xz uτ xx + vτ xy + wτ xz + q x       S(W c ) =       0 τ yx τ yy τ yz uτ yx + vτ yy + wτ yz + q y       T (W c ) =       0 τ zx τ zy τ zz uτ zx + vτ zy + wτ zz + q z      
To close the equations, the pressure is modeled by the perfect gas state law E = p (γ-1)ρ + 1 2 V • V , the heat flux q is modeled by using a Fick or Fourier law q(ǫ, κ) = -γµ Pr ∇ǫ, the adiabatic index is set to γ = 1.4, the Prandtl number is set to Pr = 0.72 and the viscosity µ is either assumed to be constant or modeled with the Sutherland law. By neglecting the viscous fluxes and the source term in the equation 1, the Euler equations are recovered.

Spatial

discretization A Mixed finite-Element/finite-Volume (MEV) discretization is used, which consists in discretizing the domain with a mixed finite-element/finite-volume approach of vertex centered type. The inviscid fluxes are discretized with a finite-volume approach while the viscous fluxes are discretized with a finite-element approach [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF][START_REF] Francescatto | M éthodes multigrilles par agglom ération dirrectionnelle pour le calcul d' écoulements turbulents[END_REF].

A polygonal bounded domain Ω ⊂ R n is considered with a bound Γ, sub-divided into a tetrahedrization or triangulation T h with elements T i . Around each vertex s i a finite-volume control cell C i of a measure m(C i ) is constructed. The set of vertices which are joined to the vertex s i is denoted by N (s i ). The subset of all the highest topological dimension polygons sharing the vertex s i is denoted by T (s i ).

The inviscid fluxes are computed on the dual control cells C i while the viscous fluxes are computed on the elements T i . A weak formulation of the Navier-Stokes equations can be expressed using a Galerkin approach. By integrating the equation 1 over a control cell S i (dual control cell or element) against a regular test function ϕ i , the weak formulation is written as:

Si (∂ t W c + ∇ • F (W c ))ϕ i dΩ = Si (∇ • N (W c ))ϕ i dΩ + Si (S(W c ))ϕ i dΩ (2)
The finite-volume method can be interpreted as a Galerking method with the control cell S i = C i and with the test function equals 1 inside the dual control cell and 0 outside. This test function, related to the control cell C i is defined as:

ϕ C i (x) = 1, if x is in C i 0, else
The variables W c are considered to be constant on each control cells C i . These constants are denoted by W c i on the cell C i (see Fig. 5 and Fig. 6).

By using the Green-Ostrogradski theorem, the left term of 2 becomes:

Ci (∂ t W c + ∇ • F (W c ))ϕ C i dΩ = m(C i )∂ t W c i + ∂Ci (F (W c ) • η))dσ (3)
with η the outward unit normal on ∂C i , the boundary of the cell C i . Furthermore, as the domain is discretised,

∂Ci (F (W c ) • η))dσ = sj ∈N (si) ∂Ci∩∂Cj (F (W c ) • η))dσ ≈ sj ∈N (si) F (W c ) |ij • ∂Ci∩∂Cj ηdσ = sj ∈N (si) F (W c ) |ij • η ij (4)
The F (W c ) |ij part is modeled by using an approximate Riemann solver for which the associated numerical flux is Φ ij = Φ(W c i , W c j , η ij ). In the current work, the Rusanov, Steger-Warming [START_REF] Steger | Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods[END_REF][START_REF] Batten | On the choice of wavespeeds for the hllc riemann solver[END_REF][START_REF] Toro | Restoration of the contact surface in the hllriemann solver[END_REF] numerical fluxes are used.

To reach a high-order approximation in space, a MUSCL reconstruction technique is used. The reconstructed primitive state from W p i at the common interface of the cells C i and C j is denoted by

W p ij : W p ij = W p i + 1 2 α ij (∇W p i ) • IJ . The slope (∇W p i ) • IJ is approximated by ∆W p i = 2/3∆ |N W p i + 1/3∆ |C W p ij .
The nodal slope ∆ |N W p i is computed from the nodal P1-Galerkin gradients that is the average gradient of the gradients computed on the elements T ∈ T (s i ). The slope ∆ |C W p ij corresponds to the centered slope W p j -W p i . We denote α ij (∆ |N W p i , ∆ |C W p ij ) the limiting coefficient. Note that the reconstruction is performed using primitive variables and not conservative variables.

Finally, the high-order numerical flux becomes:

Φ ij High order = Φ(W c ij , W c ji , η ij ).
The finite-element method can be retrieved by using the control cell S i = T i and by using basis functions on the element as test functions. The polygon T i is defined as a Lagrange P 1 finite element with a canonical basis denoted by ϕ T i for each vertex s i . The local basis function ϕ T i is equal to 1 at vertex s i and zero at the other vertices of the element. A P 1 approximation of any function f on a polygon T i can be done by projecting it on the canonical basis ϕ T of the polygon:

f (X) ≈ f Ti h (X) =< f, ϕ T >= Ns∈Ti i=1 f (i)ϕ T i (X)
A linear approximation on each polygon T i is selected by using linear polynomial functions defined by the set P 1 . This approximation is used for the density, the velocity and the temperature (or equivalently the internal energy).

By integrating by part, and by using the Green-Ostrogradski theorem, the diffusive term of 2 becomes:

Ti (∇ • N (W c ))ϕ T i dΩ = - Ti N (W c ) • ∇ϕ T i dΩ + ∂Ti (N (W c ) • n)ϕ T i dσ (5)
The velocity gradient being assumed to be constant by element, and ∇ϕ T i being a constant:

Ti N (W c ) • ∇ϕ T i dΩ = N (T i ) • ∇ϕ T i Ti dΩ = N (T i ) • ∇ϕ T i m(T i ) (6)
By denoting η Ti = -ϕ T i m(T i ), the viscous numerical flux is defined as:

Υ Ti = N (T i ) • η Ti (7)
Thus the equation 5 becomes:

Ti (∇ • N (W c ))ϕ T i dΩ = Υ Ti + ∂Ti (N (W c ) • n)ϕ T i dσ (8)
The term ∂Ti (N (W c ) • n)ϕ T i dσ concerns the boundary conditions for the viscous terms.

Time integration An implicit second-order time discretization is obtained by using a dual time step approach and a backward time integration. The dual time step approach is similar to the one proposed in [START_REF] Jameson | Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings[END_REF]. We introduce the computed residual R i for the control cell i. By denoting δ 2 1 Λ = Λ 2 -Λ 1 the variation of the variable Λ from the state 1 to the state 2, we obtain as second-order implicit time integration scheme:

m(C i ) 3W c i n+1 -4W c i n + W c i n-1 2δ n+1 n t + R i (W cn+1 ) = 0 (9)
To solve this problem we introduce a subiteration state of index k, such as:

W c i n+1 = lim k→∞ W c i (k+1)
The linearization of 9 around the state of index k, with a local time step δ k+1 k t i yields:

m(C i ) δ k+1 k t i + 3m(C i ) 2δ n+1 n t I n + J * (W c(k) ) δ k+1 k W c i = -R i (W c(k) ) + m(C i ) δ n n-1 W c i -3δ k n W c i 2δ n+1 n t (10)
with J * (W c(k) ) an approximate Jacobian of the numerical fluxes. It is composed of an inviscid and a viscous part. The inviscid Jacobian is based on the first-order Rusanov flux. The use of the Rusanov flux, or spectral radius Jacobian approximation, is usually used in matrix-free approaches [START_REF] Kloczko | D éveloppement d'une m éthode implicite sans matrice pour la simulation 2D-3D des écoulements compressibles et faiblement compressibles en maillages nonstructur és[END_REF].

The viscous Jacobian is based on the exact Jacobian of the viscous fluxes and is computed as in [START_REF] Fezoui | Resolution numerique des equations de Navier-Stokes pour un fluide compressible en maillage triangulaire[END_REF].

The resulting Jacobian matrix is inversed with either a Jacobi or a symmetric Gauss-Seidel iterative algorithm. The use of a delta form for the implicit formulation allows to use an approximated Jacobian without loosing the second order accuracy in space reached with MUSCL extrapolation.

Boundary conditions

The boundary conditions are implemented in a weak form through numerical fluxes. The representation of this implementation is given on the Fig. 7.

W c b = W c b (W c i , W c o )
. In this study, the computation of the fictive state is carried out using the Riemann invariants for free-stream conditions, imposed pressure with extrapolation of characteristics for subsonic outlets (characteristic based methods), and symmetry conditions for slipping walls.

The non-slip condition is strongly imposed by setting a zero velocity on the boundary. The Jacobian matrix is thus modified such that the boundary velocity does not evolve.

The viscous boundary fluxes are neglected for free-stream boundary conditions and computed for other cases.

Possible implementation levels

The numerical methods presented above can be implemented in the platform and plugins with different integration levels. On the one side, we can construct a simple composition based on sophisticated plugins, that basically represent the mesh generator, the solver and the view. On the other side, we can also construct a sophisticated composition based on simple plugins, that represent all single operations performed during the computational scenario. Using the former approach the visual programming tool is not much employed, whereas using the latter approach it is intensively used.

Actually, we choose an intermediate level, that allows to obtain a modular implementation where needed by our research topic. The whole mesh generator is embedded in a single plugin, whereas the solver is partitioned into elementary plugins. Therefore, the loops over finite-volume control cells or finite-elements are defined by the composer, while the local computations of fluxes, Jacobians, time-steps, etc. are embedded into elementary plugins. We obtain finally a quite sophisticated composition, presented in Fig. 8, which contains sub-compositions, such as the one illustrated in Fig. 9 used to compute and store the inviscid flux and the local time step for a given pair of control cells.

Interaction parameters

As explained above, the user can interact with the computation by modifying at runtime a plugin or its parameters. In the case of compressible flow simulations, we have so far used three different ways to modify interactively the computations:

• Modification of a plugin to change the algorithm used. For instance, we change the numerical flux or the limiting coefficient function;

• Modification of a parameter of a plugin to change boundary conditions. For instance, we change free-stream velocity, or pressure;

• Modification of a parameter of a plugin to change a numerical method. For instance, we increase or decrease the CFL number used to compute the time step.

All these modifications are performed using a graphical interface in the composer space. Illustrations are provided in the next sections.

SOME ILLUSTRATIONS

Sonic boom crossing for F16 aircraft

As first illustration, we consider the threedimensional inviscid compressible flow around a F-16 aircraft. The computational scenario chosen here, with which the user interacts, is not the resolution of the state equation, but the construction of a reduced order model from a flow database. A set of 16 computations have been performed and stored, with a grid of about 230,000 nodes, for different values of incidence and free-stream velocity. Then, the visual programming tool is used to construct a scenario that reads these solutions, interpolates linearly the solutions to compute flow variables for a given angle of attack and free-stream velocity, and visualize the data of interest. Using some widgets, the user can modify the angle of attack and the free-stream velocity, while the visualization is updated interactively. To obtain a suitable perception of 3D fields, a volumic rendering feature is implemented in the platform, which lets the user control the color and transparency associated to each field value. The scalar field of interest is evaluated on a cartesian grid from the flow solution and updated as the user modifies the interaction parameters. This task is carried out using GPUs, in order to obtain a satisfactory rendering. The pressure isolines on the aircraft surface in subsonic regime are illustrated in Fig. 10. The use of the volume rendering technique to obtain a 3D perception of the pressure field is shown on Fig. 11. As the user increases the free-stream velocity, the flow becomes supersonic and three shock waves appear, on the nose, the wing and the tail of the aircraft (see Fig. 12). It is possible to modify the lookup table used for the volumic rendering to isolate the different shock waves, as illustrated in Fig. 13.

Supersonic ramp

In the previous example, the user interacts with a reduced order model of the flow fields. Now, we would like to illustrate the possible interaction with the PDE solver itself. We consider the inviscid supersonic flow over a 15 • ramp, for inlet conditions corresponding to a Mach number of value M in = 2. The grid counts about 2,500 nodes. An unsteady simulation is performed in order to observe the flow development. During a first phase, one can observe the development of a shock wave and a rarefaction wave (see Fig. 14). Then, the shock is reflecting on the opposite wall and the flow converges to the solution depicted in Fig. 15. At this time, the user modifies the inlet boundary conditions using the widgets located at the left of the visualization space: the Mach number is progressively increased to M in = 3. The inlet flow change is progressing in the tunnel and modifies the shock wave characteristics, as illustrated in Fig. 16. This experiment illustrates how the user can interact with a computation at runtime and visualize simultaneously the flow modifications. The change of boundary condition value has been demonstrated here, but we can consider other types of interaction: the user can for instance modify the boundary condition type or location, or even modify numerical parameters, such as the numerical flux evaluation or the time step.

Oscillatory jet actuation

Finally, we present the interactive simulation of a viscous compressible flow. We consider the laminar flow over a flat plate, well known as Blasius test-case, and we introduce an oscillatory suction/blowing jet. The free-stream Mach number is M in = 0.3 and the Reynolds number is Re = 100 (based on the distance δ between the leading edge and the jet location) . The oscillatory jet width is δ/10 and frequency f = 10. The grid counts about 8,000 nodes. A second-order time integration is employed, with a physical time step 5 10 -4 . During a first phase of the simulation, no actuation is done. A boundary layer is developing, as illustrated in Fig. 17. Then, at a time t 0 , the user decides to activate the actuation. This is done by setting a non-zero actuation amplitude for the plugin that computes the jet boundary conditions (see Fig. 18). For the next time steps, the oscillatory blowing / suction boundary conditions modify significantly the boundary layer, as can be seen in Fig. 19 to Fig. 22. Although this is not presented here, the user can also modify the actuation frequency interactively.

DISCUSSION AND CONCLUSION

A modern software architecture, based on a platform / plugins system and a visual programming tool, has been used to implement numerical methods for compressible flow simulations. The platform features allow the user to modify interactively the computations by changing some parameters or even some plugins at runtime, while observing the impact on the solution.

Some tests have been carried out, that deal with reduced-order model construction for the 3D inviscid flow around an aircraft, interaction with inlet boundary conditions for the supersonic flow over a ramp, and finally modification of an oscillatory jet for a boundary layer flow.

These tests have shown that interactive computation and visualization in CFD is possible and is beneficial from a scientific point of view:

• a better understanding of physical phenomena or numerical behaviors is obtained;

• the use of CFD tool is more attractive for nonexperts;

• the adjustment of physical or numerical parameters is easier.

Of course this approach is limited by the computational facility used. Presently, for a sequential approach, the interactive computation is satisfactory until some dozen of thousand nodes. Therefore, we are presently implementing the use of parallel approaches in the platform, to be able to apply interactive computation to large-scale problems. Our objective is to study interactively turbulent flows, by using simultaneously high-performance computing and virtual reality facilities.

Figure 1 :

 1 Figure 1: Illustration of plugin system.

Figure 2 :

 2 Figure 2: Illustration of the visual programming system.

Figure 3 :

 3 Figure 3: Illustration of the visualization space.

Figure 4 :

 4 Figure 4: Virtual reality facility at INRIA Sophia Antipolis M éditerran ée center.

Figure 5 :

 5 Figure 5: Illustration of a control cell in 2D.

Figure 6 :

 6 Figure 6: Illustration of a control cell in 3D.

Figure 7 :

 7 Figure 7: Illustration of the implementation of weak boundary conditions.

Figure 8 :

 8 Figure 8: Composition used for flow simulation.

Figure 9 :

 9 Figure 9: Sub-composition used for flux and local time step computation and storage.

Figure 10 :

 10 Figure 10: Pressure isolines for F16 aircraft.

Figure 11 :

 11 Figure 11: Pressure field in subsonic regime.

Figure 12 :

 12 Figure 12: Pressure field in supersonic regime.

Figure 13 :

 13 Figure 13: Modification of the volume rendering lookup table to isolate shock waves.

Figure 14 :

 14 Figure 14: Pressure field: shock wave and rarefaction wave in development.

Figure 15 :

 15 Figure 15: Pressure field: shock wave reflecting.

Figure 16 :

 16 Figure 16: Pressure field: interaction between inlet condition change and shock wave.

Figure 17 :

 17 Figure 17: Velocity modulus without actuation (t 0 ).

Figure 18 :

 18 Figure 18: Velocity modulus (t 0 + 0.015).

Figure 19 :

 19 Figure 19: Velocity modulus (t 0 + 0.040).

Figure 20 :

 20 Figure 20: Velocity modulus (t 0 + 0.060).

Figure 21 :

 21 Figure 21: Velocity modulus (t 0 + 0.085).

Figure 22 :

 22 Figure 22: Velocity modulus (t 0 + 0.120).


Louis Blanchard

Régis Duvigneau

Bernard Mourrain

Automated Generation of Geometric Models for a Bulbous Trawler Ship

Keywords: parameterization, shape, bow, splines * Opale Project-Team paramétrisation, forme, bulbe, splines

The objective of this work is to develop an automated geometric modeler for trawler ships including a bow, for shape optimization purpose. Two approaches are tested to generate a spline-based representation of the hull. For the first one, a starting shape without bow is deformed to reach a target composed of points or lines, while preserving the shape regularity using a penalty functional. For the second approach, a starting hull including a baseline bow is deformed by handling meta-parameters, such as bow length, thickness or orientation. The different approaches are finally illustrated for a realistic application.

Introduction

Numerical simulations associated with automated design optimization tools can now be used for practical problems in naval hydrodynamics, thanks to the robustness of simulation methods and the increase of computational facilities. A major difficulty in this context is the automated generation of geometric models, defined by a few design parameters, and the corresponding computational grid.

Therefore, we study in this work how to generate automatically a parameterized hull shape model, including a bow, with application to trawler ships. The final objective is to include the developed tool into a design optimization procedure, that allows to automatically optimize the bow shape.

The primary objective of the geometrical model is to facilitate the modification of an existing hull shape by modifying a small number of parameters, that are meaningful for naval architects (see figure 1).

Two approaches are considered successively: for the first one, a starting shape without bow is deformed to reach a target composed of points or lines, while preserving the shape regularity using a penalty functional. For the second approach, a starting hull including a baseline bow is deformed by handling meta-parameters, such as bow length, thickness or orientation. 2 Fitting approach

Distance error minimization

We propose first to generate bow shapes from an existing hull without bow. The method consists in fitting the hull surface to a target bow defined by point clouds, by distance minimization.

We refer to articles [START_REF] Chen | Computing the minimum distance between a point and a clamped B-spline surface[END_REF], [START_REF] Helmut Pottmann | Approximation with active bspline curves and surfaces[END_REF] and [START_REF] Wang | Fitting b-spline curves to point clouds by squared distance minimization[END_REF] for a detailed review of distance minimization problems. Three methods are commonly used to approximate a target shape defined by a set of points cloud: the point distance minimisation (PDM), tangent distance minimisation (TDM) and Squared Distance minimization (SDM). See article [START_REF] Wang | Fitting b-spline curves to point clouds by squared distance minimization[END_REF] and the internship report [START_REF] Proust | Modelisation et optimisation de carenes[END_REF] for a comparison of these methods.

Definition of distance error term

We denote ⊗ the tensorial product. Let us consider for ξ ∈ R 3 the B-spline surface:

σ c (ξ) def = Id 3 ⊗ B T (ξ) c (1)
where B(ξ) ∈ R n is the vector of B-spline basis functions and c ∈ R 3n is the vector of control points.

We consider the following problem: for a given bow shape, defined by a set of data points X k , ∀k ∈ [1, M ], compute the B-spline surface to approximate the points X k .

For a fixed k ∈ [1, M ] we denote σ c (ξ k ) the point of the surface σ c (ξ) that minimizes the distance between the point X k and the surface σ c (ξ):

d σ c (ξ k ), X k = argmin σ c (ξ) -X k R 3 . Using the following notations : B k def = B(ξ k ) , A k def = Id 3 ⊗ B T k we get: σ c (ξ k ) = Id 3 ⊗ B T k c = A k c (2)
We defined the Point Distance minimization (PDM) function as follows :

e k (c) = 1 2 σ c (ξ k ) -X k 2 R 3 (3)
We defined the Tangente Distance minimization (TDM) function as follows :

ẽk (c) = 1 2 σ c (ξ k ) -X k , N k R 3 2 (4)
where N k is the external normal vector of the point σ c (ξ k ).

Finally, we defined the Squared Distance minimization (SDM) function as follows :

ēk (c) = ẽk (c) + 1 2 d k d k +ρ u k σ c (ξ k ) -X k , T u k R 3 2 + d k d k +ρ v k σ c (ξ k ) -X k , T v k R 3 2 if d k > max(ρ u k , ρ v k ) 0 otherwise (

Minimization

In this section we consider the minimization of the distance error term. The problem can be seen as the non linear problem:

J(c) = M k=1 e k (c) , J(c) = M k=1 ẽk (c) , J(c) = M k=1 ēk (c) (6)
To solve these problems, we consider deterministic descent methods, based on the computation of gradient and of Hessian of the cost function. The iterative method consists in computing the update c t+1 = c t + α t d t , where d t ∈ R n is a descent direction and α t ∈ R + is a step length obtained by a line search :

α t = arg min α∈R + J(c t + α d t ) (7)
We refer to the book [START_REF] Gill | Practical optimization[END_REF] for a detailed review of deterministic optimization methods. Here, we simply consider the following methods :

• gradient method : with the following update for the descent direction :

∀t ∈ N * , d t = -∇J(c t ) (8)
• Newton-Raphson method : where the descent direction is the solution of the linear system :

∀t ∈ N , ∇ 2 J(c t ) d t = -∇J(c t ) (9)
The next propositions allow to evaluate the derivatives of cost functions defined by Eq [START_REF] Wang | Fitting b-spline curves to point clouds by squared distance minimization[END_REF].

Proposition: The fist and second derivatives of the cost function Eq(6) based on the Point Distance minimization function Eq(3) are :

∇J(c) = M k=1 ∇e k (c) where ∇ e k (c) = σ c (ξ k ) -X k ⊗ B k ∇ 2 J(c) = M k=1 ∇ 2 e k (c) where ∇ 2 e k (c) = Id 3 ⊗ B k B T k (10)
Proof: Using the definition of Point Distance minimization function from Eq(3), we get

e k (c) = 1 2 A k c -X k , A k c -X k R 3 = 1 2 c , A T k A k c R 3n -2 c , A T k X k R 3n + X k 2 R 3
Consequently, the gradient of the Point Distance minimization function turns into :

∇ e k (c) = A T k A k c -A T k X k = Id 3 ⊗ B k Id 3 ⊗ B T k c -Id 3 ⊗ B k X k = Id 3 ⊗ B k B T k c -Id 3 ⊗ B k X k = Id 3 ⊗ B k σ c (ξ k ) -Id 3 ⊗ B k X k = σ c (ξ k ) -X k ⊗ B k
and the Hessian of the Point Distance minimization function turns into :

∇ 2 e k (c) = A T k A k = Id 3 ⊗ B k Id 3 ⊗ B T k = Id 3 ⊗ B k B T k RR n°8011
Proposition: The fist and second derivatives of the cost function Eq(6) based on the Tangente Distance minimization function Eq(4) are :

∇ J(c) = M k=1 ∇ẽ k (c) where ∇ ẽk (c) = σ c (ξ k ) -X k , N k R 3 N k ⊗ B k ∇ 2 J(c) = M k=1 ∇ 2 ẽk (c) where ∇ 2 ẽk (c) = N k N T k ⊗ B k B T k (11)
Proof: Using the definition of the Tangente Distance minimization function from Eq(4), we get

ẽk (c) = 1 2 c , A T k N k R 3n -X k , N k R 3 2
Consequently, the gradient of the Tangente Distance minimization function turns into :

∇ ẽk (c) = σ c (ξ k ) -X k , N k R 3 A T k N k = σ c (ξ k ) -X k , N k R 3 Id 3 ⊗ B k N k = σ c (ξ k ) -X k , N k R 3 N k ⊗ B k
and the Hessian of the Tangente Distance minimization function turns into :

∇ 2 ẽk (c) = (A T k N k ) (A T k N k ) T = A T k N k N T k A k = Id 3 ⊗ B k N k N T k Id 3 ⊗ B T k = N k ⊗ B k N T k ⊗ B T k = N k N T k ⊗ B k B T k Proposition:
The fist and second derivatives of the cost function Eq(6) based on the Squared Distance minimization function Eq(5) are :

∇ J(c) = M k=1 ∇ē k (c) , ∇ 2 J(c) = M k=1 ∇ 2 ēk (c) (12)
where

∇ ēk (c) = ∇ ẽk (c) + d k d k +ρ u k σ c (ξ k ) -X k , T u k R 3 T u k + d k d k +ρ v k σ c (ξ k ) -X k , T v k R 3 T v k ⊗ B k 0 (13) ∇ 2 ēk (c) = ∇ 2 ẽk (c) + d k d k +ρ u k T u k T u k T + d k d k +ρ v k T v k T v k T if d k > max(ρ u k , ρ v k ) 0 otherwise ( 14
)
Proof: Using that d k is supposed to be fixed on c, and proposition Eq(10)

Numerical experiments

We present in this section the results obtained using the Newton-Raphson method from equation Eq(9). The cost function considered is J(c) defined by Eq [START_REF] Wang | Fitting b-spline curves to point clouds by squared distance minimization[END_REF]. The initial hull shape, without bow, originates from the test-case defined in the national project BULBE (see acknowledgement Inria section). The bow shape, defined by point clouds X k ∀k ∈ [1, M ] is represented by two line targets (in red on figure 2). Figure 2 represents, initial shapes of the trawler ship (on the left) and optimal shapes of the trawler ship (on the right). These tests illustrate the efficiency of the Distance error term minimization approach applied to a trawler ship for generating bows (see figure 2). However, we have observed that this approach suffers from serious drawbacks: the method does not care about the smoothness of the resulting shapes. In practice, the bow shapes obtained exhibit irregularities that are not suitable for industrial use. To overcome this difficulty, we propose to introduce a penalty term that will regularize the problem and ensure a smooth solution surface.

Smoothness measure minimization

Definition of smoothness measure term

The smoothness measure term of a B-spline surface S = σ c (ξ) is the space integral over the B-spline surface S of the first order derivatives of the mapping function:

J s (c) = σc(ξ) tr Dσ c (ξ)[Dσ c (ξ)] T dΩ = σc(ξ) ∂ u σ c (ξ) 2 R 3 + ∂ v σ c (ξ) 2 R 3 dΩ (15)

Optimization of smoothness measure term

Proposition: The derivatives of smoothness measure from Eq(15) are:

∇ J s (c) = S ∂ u σ c (ξ) ⊗ B u (ξ) + ∂ v σ c (ξ) ⊗ B v (ξ) dΩ ∇ 2 J s (c) = Id 3 ⊗ S B u (ξ)[B u (ξ)] T + B v (ξ)[B v (ξ)] T dΩ (16)
where :

B u (ξ) = ∂ u B(ξ) , B v (ξ) = ∂ v B(ξ) (17)
Proof: Using equations Eq(1), Eq(17) and the following noation

A u (ξ) = ∂ u A(ξ) and A v (ξ) = ∂ v A(ξ), we get ∂ u σ c (ξ) = Id 3 ⊗ [B u (ξ)] T c = A u (ξ) c ∂ v σ c (ξ) = Id 3 ⊗ [B v (ξ)] T c = A v (ξ) c (18)
Consequently the definition of the smoothness measure from Eq(15) becomes:

J s (c) = S c , [A u (ξ)] T A u (ξ) c R 3n + c , [A v (ξ)] T A v (ξ) c R 3n dΩ = c , S [A u (ξ)] T A u (ξ) + [A v (ξ)] T A v (ξ) dΩ c R 3n
Consequently, the gradient of the smoothness measure from Eq(15) turns into :

∇ J s (c) = S [A u (ξ)] T A u (ξ) + [A v (ξ)] T A v (ξ) dΩ c = S [A u (ξ)] T ∂ u σ c (ξ) + [A v (ξ)] T ∂ v σ c (ξ) dΩ = S Id 3 ⊗ B u (ξ) ∂ u σ c (ξ) + Id 3 ⊗ B v (ξ) ∂ v σ c (ξ) dΩ = S ∂ u σ c (ξ) ⊗ B u (ξ) + ∂ v σ c (ξ) ⊗ B v (ξ) dΩ
and the hessian of the smoothness measure from Eq(15) turns into :

∇ 2 J s (c) = S [A u (ξ)] T A u (ξ) + [A v (ξ)] T A v (ξ) dΩ = S Id 3 ⊗ B u (ξ)[B u (ξ)] T + B v (ξ)[B v (ξ)] T dΩ = Id 3 ⊗ S B u (ξ)[B u (ξ)] T + B v (ξ)[B v (ξ)] T dΩ Inria 2.

Numerical experiments

This section presents some numerical experiments using an abstract three dimensional surface, to demonstrate the efficiency of the Newton-Raphson optimization method, applied to the reduction of the smoothness measure term J s from equation Eq(16). The initial shape exhibits very large irregularities, that we would like to reduce. It is represented in blue on figure 3 and is composed of 20 2 control points. The boundaries points define a square and are supposed to be fixed during the optimization process. Consequently the number of control points to be optimized is 18 2 (in three dimension). Consequently the size of the optimization problem is 972.

The figure 3 illstrates the evolution of the surface, for different iterations of the optimization process. The figure 4 presents the evolution of the cost function J s (in blue) and the norm of the cost function |∇J s | (in red) during the optimization process.

Figure 3: Initial shape (blue) shape after 400 iterations (red), shape after 600 iterations (yellow) and optimal shape (grey).

We obtain promising results for the Smoothness measure term minimization approach applied to a three dimensional surface test-case (see figures 3 and4).

However, for the initial problem which consists in generating a bow shape, we have to consider an objective function that combines both the Distance error term and the Smoothness mesure term:

J(c) + λ J s (c) ( 19
)
where λ is a positive coefficient used to modulate the weight of J s . It has been found in practical tests with the trawler ship that the results depend strongly on the penalty coefficient and the application of this approach is tedious for industrial cases. In this section we propose to generate bow shapes by modifying some meta-parameters, such as width, length or leading direction, from a baseline bow shape.

Using the geometrical model of the trawler ship without bow, we extract a sub-surface that would correspond to the location of a bow surface (see figure 5 on the left). Then, we generate a baseline bow surface by moving the control points of this sub-surface (see figure 5 on the right). The whole geometrical model (hull+bow) can finally be obtained by reconnecting the two surfaces.

The geometrical characteristics of the bow can be controlled by modifying some metaparameters, such as width, length or leading direction. Each modification of a parameter value is applied to the control points that define the bow surface, as a global transformation. Figure 6 presents some different bow shapes obtained by modifications of width and thickness parameters.

Obviously, this approach is far easier to use than the previous one, but it is more limited in terms of resulting shapes. For design optimization purpose, the geometric modeler should be embedded into the design procedure which is composed of the grid generation tool, the Computation Fluid Dynamics (CFD) solver and the optimization tool.

In the framework of the national project BULBE , we intend to use the commercial software package FINE/Marine from NUMECA, that gathers the grid generation software Hexpress and the solver Isis-CFD, in collaboration with K-Epsilon company.

The geometrical model, based on the methods presented above, has been implemented in the AXEL Platform developed by INRIA Galaad Project-Team. To export the resulting geometry to the grid generation tool, some additional steps have been introduced: firstly, the hull surface is embedded into a closed box used to define the computational domain. Secondly, all the surfaces are triangulated to obtain a discrete description of the domain boundary, written in an input file for the grid generation tool. Figure 7 illustrates these two steps necessary for exporting the geometrical model. We refer to [START_REF] Perrier | Module de maillage parametrique pour la cfd[END_REF] (C annex) for details concerning the description of the file.

Conclusion

Two approaches have been developed and tested to generate automatically bow shapes for a trawler ship.

The first one, based on fitting methods with regularization, seems to be promising because it allows to use very detailed shape targets. However, this approach still requires to mature for industrial use. In particular, one should improve the constraints applied on the shape in order to avoid unrealistic oscillations.

The second one, based on a direct handling of meta-parameters, is far easier to use but is more limited in terms of shape definition.

Figure 1 :

 1 Figure 1: Initial shape (left) and deformed shape to generate a bow (right) for a trawler ship.

2 R 3

 23 5) where T u k and T v k are the tangential vectors of the point σ c (ξ k ), (ρ u k , ρ v k ) are the curvatures of the point σ c (ξ k ) and where d k = e k (c) = σ c (ξ k ) -X k is supposed to be fixed on c.

Inria

Figure 2 :

 2 Figure 2: Initial shape (left) and deformed shape to generate a bow (right) for the trawler ship, to reach two line targets (in red).

Figure 4 :

 4 Figure 4: Evolution of the cost function during the optimization process.

Figure 5 :

 5 Figure 5: Initial shape (left) and deformed shape to generate a bow (right) for the trawler ship.

Figure 6 :

 6 Figure 6: Initial shape (top), deformed shapes : length = 1.5 width = 0.5 (middle), length = 2.5 width = 0.8 (bottom) .

Figure 7 :

 7 Figure 7: Domain description for computation (the shape of the trawler ship and the box) .

RR n°8011

Acknowlegement

This work is funded by the Ministry of Fishing DPMA (Direction des Pêches Maritimes et de l'Aquaculture) and the European Fund for Fishing. It gathers INRIA OPALE Project-Team, K-Epsilon company (specialized in CFD for naval hydrodynamics, see http://www.k-epsilon. com) and PROFIL company (naval architect, see http://www.jean-et-frasca.com). The objective is to design and optimize bow shapes for trawler ships, in order to reduce the fuel consumption during fishing campaigns.

Inria


J N Mei

Jean-Louis L Soubeyroux

Jean-Jacques J Blandin

J S Li

H C Kou

H Z Fu

Lian Zhou

Structural relaxation of Ti40Zr25Ni8Cu9Be18 bulk metallic glass

Keywords: Bulk metallic glass, Titanium alloy, Structural relaxation, Quenched-in nuclei, Glass forming ability, Thermal stability

HAL is

Introduction

Since a glass can be described as an undercooled melt, such glassy structure is thermodynamically metastable compared with the ideal glassy state and tends to undergo structural rearrangement to reach the ideal glassy state (structural relaxation) [START_REF] Chen | Glassy metals[END_REF]2] or transform to an equilibrium crystalline structure through a series of continuous transformations once the external energies are supplied either thermally by heat treatments or mechanically by application of external stresses, leading to continuous changes in physical, chemical and mechanical properties.

Ti-based bulk metallic glasses (BMGs), as promising lightweight structural materials, have attracted more and more attention during recent years [3][4][5][6][7][8][9][10][11][12]. Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG, as one of the best Ti-based BMGs known today, was found to exhibit a good combination of glass forming ability (GFA, the maximum diameter could be 8 mm) and mechanical properties (high compressive strength and large compressive plastic strain) [9,13,14], which could be a potential structural material. Our recent study has investigated its crystallization behavior, and found its unique quenched-in nuclei/amorphous matrix structure [14,15], suggesting that this alloy is a marginal BMG. Actually such marginal glassy alloys exist extensively and are complex with respect to typical fully glassy alloys, accordingly may lead to the different thermally induced transformation behaviors for their different structural nature and the microstructure sensitivity to thermal activation. However, so far there is no detailed investigation on structural relaxation of such complex marginal BMGs. Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG could be a good medium for investigating the structural relaxation behavior of such complex metallic glasses. Therefore, the present work aims to investigate the structural relaxation behavior of Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG.

Experiment

The ingots of master alloy with nominal composition of Ti 40 Zr 25 Ni 8 Cu 9 Be 18 (at.%) were prepared by melting elements of 99.8 to 99.999 wt.% purity with high frequency electromagnetic induction in a water-cooled copper crucible under a high purity argon atmosphere. The ingot alloys were remelted at high temperature several times to ensure the compositional homogeneity. Cylindrical samples with 2 mm in diameter and 40 mm in length were prepared by injection casting into a water-cooled copper mould via an electromagnetic levitation-melt state under a high purity argon atmosphere before casting to obtain a completely melted state. Alloy's oxygen content was measured <400 ppm.

The glass transition temperature (T g ) and the onset crystallization temperature (T x1 ) were determined to be 313 and 372 o C, respectively, by differential scanning calorimetry (DSC, Netzsch DSC 404s) measurement at a heating rate of 20 K/min in a flowing argon atmosphere. The microstructure was observed by high-resolution transmission electron microscopy (HRTEM, JEM-3010).

Structural relaxation was investigated through measuring the relaxation of the frozen-in enthalpy towards the equilibrium state as a function of time at a certain constant temperature by using DSC. The enthalpy reduction during an isothermal relaxation below the calorimetric glass transition is recovered as an "overshooting" at the glass transition during a subsequence DSC run. The detailed temperature program can be summarized as followed: to erase the thermal history of metallic glasses, each sample was heated from room temperature to a certain temperature (about 20 °C above T g ) at a heating rate of 20 K/min, and then immediately cooled down to room temperature at a cooling rate of 50 K/min, then heated the sample to the isothermal annealing temperature T (273, 283, 293 and 303 °C) at a heating rate of 30 K/min and annealed for different period of time t (5 min-6 h). After each isothermal annealing, the sample was heated up to high temperature to complete the crystallization (600 °C) at a heating rate of 20 K/min and then immediately cooled down to room temperature at a heating rate of 50 K/min for the following calculation of recovery enthalpy. Finally, this heating and cooling cycle was repeated using the same procedure as before for baseline correction. The relaxation recovery enthalpy (ΔH T (t)) was determined by measuring the overshoot endothermic enthalpy around T g of an annealed sample with respect to the as-cast sample. It is well known that the time dependence of ΔH T (t) at a given temperature of amorphous materials generally can be fitted by a stretched exponential relaxation function [16][17][18] ]}

Results

) ( exp[ 1 { ) ( β τ t H t H eq T T - - Δ = Δ ( 1
)
where τ is the mean structural relaxation time, eq T H Δ

is the equilibrium value of ΔH T (t) for t→∞ at the annealing temperature T, β is the stretching exponent whose value is between zero and unity 0<β≤1, reflecting the width of relaxation time distribution spectrum originating from the inhomogeneity of the supercooled liquid. The value of β=1 corresponds to a single relaxation time with exponential behavior.

The smaller the value of β, the more the distribution of atomic motions deviates from a single exponential behavior. So, if the relaxation time distribution becomes broader, the value of β will be reduced, and the enthalpy recovery curve will be spread out over a wider temperature range. ) exp(

0 0 0 T T T D - = * τ τ (2)
where τ 0 is the pre-exponential factor related to relaxation time constant, D * is the fragility parameter and T 0 is the VFT temperature at which barriers with respect to flow would go to infinity and is believed to correspond to the theoretical Kauzmann temperature T k . In most enthalpy relaxation studies, T 0 in VFT equation is assumed to be 0 and accordingly the mean relaxation time τ for the enthalpy relaxation is assumed to follow an Arrhenius behavior [21] ) exp(

' 0 a a RT E τ τ = (3)
where ' 0 τ is the pre-exponential factor, R is the ideal gas constant, E a is the apparent activation energy [18,20]. It is found that the experimental annealing temperature dependence of the mean relaxation time τ for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG can be well fitted with an Arrhenius plot as shown in Fig. 2, and yields a ' 0 τ of 1.96×10 -19 s and an E a of 233 kJ/mol (2.42 eV). It should be pointed out that such Arrhenius behavior cannot hold over a wide temperature since that actually T 0 is not equal to 0 and approximately several tens of Kelvins below T g for glasses. However, the T 0 item in VFT equation can be ignored when considering the temperature dependence of τ in the experimental annealing temperature range.

Discussion

Although the unique quenched-in nuclei/amorphous matrix structure, surprisingly, its enthalpy relaxation behavior follows the general characteristics of amorphous materials without being disturbed The characteristic parameters, stretching exponent β and mean structural relaxation τ derived from the stretched exponential relaxation function, are essential for the analysis of enthalpy relaxation, which are believed to be correlated with the GFA and thermal stability of metallic glasses, respectively.

However, their characteristic of temperature dependence brings about difficulties for constructing the correlation between stretching exponent β and GFA, and correlation between mean structural relaxation τ and thermal stability. In this study, two parameters, β g and τ g , representing the stretching exponent and the mean structural relaxation time at the calorimetric glass transition temperature, respectively, were proposed to be the indicators of enthalpy relaxation properties.

The relation between β g and GFA of metallic glasses was constructed through fragility index m as a chain, since the stretching exponent β has a correlation with the fragility index m which is an indicator of the fragility of supercooled liquid and the strong liquid (with small fragility index) is generally correlated with the good GFA. It is known that most of the bulk metallic glasses show strong liquid behavior and generally have a fragility index of less than 50 [19,22,23]. And a reduced m=42/β relation can be derived from a relation between the fragility index m and stretching exponent β based on the defect diffusion model proposed by Bendler et al. [24], this expression has been confirmed to be suitable for describing the m-β relation for polymeric, alcoholic and metallic glasses. Accordingly, the value of stretching exponent β for the metallic glasses who exhibit a good GFA empirically should be larger than 0.84. Fig. 4 shows the stretching exponent β as a function of annealing temperature T for

Conclusions

Enthalpy relaxation from the amorphous state into the supercooled liquid state was investigated on the Ti 40 Zr 25 Ni 8 Cu 9 Be 18 bulk metallic glass below the calorimetric glass transition. The enthalpy relaxation behavior of this alloy follows the general characteristics of enthalpy relaxation for glasses although its unique quenched-in nuclei/amorphous matrix structure: a stretched exponential function with the relaxation time obeying an Arrhenius law, which is due to the incubation time for crystallization τ inc is longer than the relaxation time τ and the enthalpy relaxation is almost not disturbed by the crystallization of quenched-in nuclei. Two parameters, β g and τ g , representing the stretching exponent and the mean structural relaxation time at the calorimetric glass transition temperature, respectively, were proposed to correlate with GFA and thermal stability, respectively. The high value of β g that is much higher than 0.84 and approaches unity for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG confirms its good GFA as typical Zr-based BMGs; however, its low value of τ g (105 s) indicates a worse thermal stability compared with typical BMGs.

Fig. 1 Fig. 1

 11 Fig. 1 shows the recovery enthalpy ΔH T (t) as a function of annealing time t at various annealing temperatures T of 273, 283, 293 and 303 o C for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG. It can be seen that the value of ΔH T (t) increases rapidly with increase in t at the initial stage of relaxation, then changes to increase gradually, finally tends to saturate. Further investigation find that the ΔH T (t)-t relation does not follow a simple exponential relaxation function, which could be due to the microstructural heterogeneities of the sample, and suggests that the enthalpy relaxation of the Ti 40 Zr 25 Ni 8 Cu 9 Be 18 liquid does not involve a single Debye relaxation event, as is in accordance with the previous report related to amorphous materials.

(Fig. 3

 3 Fig. 3 Bright field image and corresponding SADP obtained from annealed Ti 40 Zr 25 Ni 8 Cu 9 Be 18 samples

Fig. 4 Fig. 5

 45 Fig. 4 Stretching exponent β as a function of annealing temperature T for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 as well as Vit 4 and Vit 106a BMGs.The relation between τ g and thermal stability of metallic glasses was constructed based on the simple idea that the long structural relaxation time corresponds to a good thermal stability. There exists a linear relationship between the logarithmic relaxation time lnτ and the difference between the glass transition temperature and annealing temperature ΔT,

Table 1 .

 1 Obviously, the experimental ΔH T (t)-t relation for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG can be well fitted with stretched exponential relaxation function. The obtained It is found that the values of Zr 46.75 Ti 8.25 Cu 7.5 Ni 10 Be 27.5 [19], Pd 43 Cu 27 Ni 10 P 20 [18] and Zr 58.5 Cu 15.6 Ni 12.8 Al 10.3 Nb 2.8[20] BMGs, indicating that the structural relaxation behavior of BMGs are complex.

	H Δ	eq T
		H Δ	eq T	, τ and β are dependent on the annealing

, τ and β, at various annealing temperatures are listed in , τ and β is mainly ascribed to the degree of the amorphous phase deviating from its equilibrium state at each annealing temperature. It should be pointed out that for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG, stretching exponent β is not constant but slightly increasing with annealing temperature. Such annealing temperature dependence of stretching exponent β has been found in some

Table 1

 1 The values of , τ and β obtained from the fitting of ΔH T (t) vs. t curves based on the stretched exponential relaxation function for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG.

	H Δ	eq T
	T (°C)	H Δ	eq T	(J/g)	τ (s)	β
	273		4.2		3997	0.86201
	283		3.7		1616	0.87987
	293		2.6		624	0.91135
	303		1.9		266	0.93578

In the literature, Vogel-Fulcher-Tammann (VFT) equation is used to describe the non-Arrhenius behavior of the mean relaxation time τ

[21]

  Fig. 2 Mean enthalpy relaxation time as a function of inverse temperature for Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG fitted with the Arrhenius equation.

	τ (s)	10 3			τ ' 0 =1.96E-19 s
					E a =233 kJ/mol
						=2.42 eV
		1.72	1.74	1.76	1.78	1.80	1.82	1.84
				1000/T (1/K)

  Ti 40 Zr 25 Ni 8 Cu 9 Be 18 BMG as well as typical Zr-based BMGs (Vit 4 and Vit 106a alloys), the β g values can be obtained through linear extrapolation. It is found that the value of β g for Ti 40 Zr 25 Ni 8 Cu 9 Be 18

	1.00
	0.95
	0.90
	0.85	Ti 40 Zr 25 Ni 8 Cu 9 Be 18		Vit 4		Vit 106a
	0.80
	0.75
	0.70	560	580	600	620	640	660	680
	BMG is comparable with other typical BMGs, much higher than 0.84 and approaches unity, implying
	its strong liquid behavior and accordingly good GFA.

Acknowledgements

This work was supported by National Natural Science Foundation of China (50601011), National Basic Research Program of China (2007CB607603), Fundamental Research Project of National Defense of China (A2720060295) and China Scholarship Council. We are very grateful to Dr. P. Wen who works in Institute of Physics (CAS) for helpful advice.