Santiago Bragagnolo
email: santiago.bragagnolo@inria.fr

Matteo Marra
email: mmarra@vub.be

Guillermo Polito
email: guillermopolito@gmail.com

Elisa Gonzalez Boix

Elisa Gonzalez
email: egonzale@vub.be

Boix Towards Scalable

Towards Scalable Blockchain Analysis

Keywords: Computing methodologies → MapReduce algorithms, • Information systems → Query languages, Information retrieval query processing, Computing platforms, Digital cash blockchain, smart contracts, big data, Map/Reduce

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

INTRODUCTION

Many existing work on blockchain analysis (on the domains of fraud, security, statistic analysis, etc.), base their algorithms and analysis strategies on analysing transactions [START_REF] Lischke | Analyzing the bitcoin network: The first four years[END_REF][START_REF] Meiklejohn | A fistful of bitcoins: Characterizing payments among men with no names[END_REF][START_REF] Vasek | There is no free lunch, even using Bitcoin: Tracking the popularity and profits of virtual currency scams[END_REF]. For gathering transaction information and running their algorithms, they gather data by sequentially selecting and filtering information from different blockchain networks. This kind of usage is becoming outdated nowadays. Ethereum and other blockchain platforms store a massive amount of heterogeneous data: transactions, accounts, smart contracts, etc. In early 2017, Ethereum was estimated to have approximately 300GB of data [START_REF] Bartoletti | A General Framework for Blockchain Analytics[END_REF], and it keeps growing at a high rate. Extracting information from this massive amount of data is not an easy task. Currently, there are two options to search for specific information stored in a block (e.g., account's data). Either the unique identifier of the block (or block hash) containing such information is known, or a sequential search in the blockchain is required starting from a block (e.g. the most recent one) and searching for every parent block for any related information. Hashes of the blocks are typically unknown (unless a secondary database is used to store all hashes of blocks deployed on the blockchain). As a result, a full scan of the blockchain may be required to be able to locate and analyse data.

In addition, Ethereum stores information (e.g. contracts, transactions) with a generic representation which does not include metadata describing the stored information. While this reduces the data size improving performance, it complicates the research of information once it is deployed in the blockchain.

The goal of our work is to ease the extraction and analysis of data from blockchain platforms like Ethereum. This is crucial to assist the process of debugging and fixing bugs in blockchain applications, particularly in smart contract applications. In prior work, we have proposed a SQL-like query language for Ethereum to ease the task of querying data in the blockchain [START_REF] Bragagnolo | Ethereum Query Language[END_REF]. Such a query language allows developers to quickly scan the blockchain to access individual elements and aggregate arbitrary values. To perform fast searches by semantic attributes, the usage of indexes is crucial. However, the setup and maintenance of indexes on the data stored in Ethereum is a highly consuming task. Some work has explored the performance problems to do blockchain analytics [START_REF] Bartoletti | A General Framework for Blockchain Analytics[END_REF][START_REF] Kalodner | BlockSci: Design and applications of a blockchain analysis platform[END_REF] but either they work on early stages of the main cryptocurrency networks (when the amount of data was still bearable) [START_REF] Bartoletti | A General Framework for Blockchain Analytics[END_REF], or on smaller blockchains than Ethereum [START_REF] Kalodner | BlockSci: Design and applications of a blockchain analysis platform[END_REF].

In this paper, we explore the usage of parallelizing techniques present in Big Data platforms, in particular Map/Reduce, to extract and analyse information from the blockchain. We focus on applying the Map/Reduce model to build indexes of blocks, transactions and contracts. We show that using our approach significantly improves the index generation time by 7.77x with a setup of 20 worker nodes running on a cluster of 10 Intel Xeon CPU E3-1240 @ 3.50GHz machines. These good results have been achieved with just a single blockchain data node running in the cluster.

In summary, our contributions are:

• to the best of our knowledge, we are the first using Big Data techniques to enable blockchain analytics. • we present a Map/Reduce architecture that can effectively and efficiently query the blockchain. • we identify the bottlenecks of such architecture and propose possible solutions for it.

STATE OF THE ART AND MOTIVATION

Prior work has identified the need of blockchain analysis for several use-cases such as the detection of attacks and security vulnerabilities [START_REF] Baqer | Stressing out: Bitcoin stress testing[END_REF][START_REF] Fleder | Bitcoin transaction graph analysis[END_REF][START_REF] Meiklejohn | A fistful of bitcoins: characterizing payments among men with no names[END_REF][START_REF] Meiklejohn | A fistful of bitcoins: Characterizing payments among men with no names[END_REF][START_REF] Ober | Structure and anonymity of the bitcoin transaction graph[END_REF][START_REF] Spagnuolo | Bitiodine: Extracting intelligence from the bitcoin network[END_REF][START_REF] Vasek | There is no free lunch, even using Bitcoin: Tracking the popularity and profits of virtual currency scams[END_REF], finantial frauds [START_REF] Moser | An inquiry into money laundering tools in the Bitcoin ecosystem[END_REF][START_REF] Möser | Towards risk scoring of Bitcoin transactions[END_REF][START_REF] Vasek | There is no free lunch, even using Bitcoin: Tracking the popularity and profits of virtual currency scams[END_REF] or just retrieving statistics or economic indicators [START_REF] Lischke | Analyzing the bitcoin network: The first four years[END_REF][START_REF] Ron | Quantitative analysis of the full bitcoin transaction graph[END_REF]. Most of these works create their custom representation of the blockchain with all the data required for the desired analysis.

To cope with blockchain analysis in a more general way, several platforms and frameworks have been proposed in literature [START_REF] Bartoletti | A General Framework for Blockchain Analytics[END_REF][START_REF] Bragagnolo | Ethereum Query Language[END_REF][START_REF] Kalodner | BlockSci: Design and applications of a blockchain analysis platform[END_REF]. Kalodner et al. [START_REF] Kalodner | BlockSci: Design and applications of a blockchain analysis platform[END_REF] propose BlockSci, a BlockChain analysis platform to analyze crypto currency based Blockchains like Bitcoin, Litecoin and others. BlockSci's main trait is its architecture: it imports the entire BlockChain data into memory, where queries can be executed fast. This makes BlockSci outperform other tooks by 15x-600x, and spends only around 60 seconds to setup their memory database because they access directly the blockchain data of big blockchains by directly accessing the raw data in disk instead of using the APIs designed for it. BlockSci, however, suffers two main drawbacks. First, BlockSci was not designed to work on smart contracts. The disregard for the smart contracts takes out their possibilities the contract based analysis. Many of the emerging blockchain technologies provide smart contract support (Ethereum1 , Hyperledger Fabric2 , NEM3 , Waves4 , etc). Second, it was not designed to work over a large amount of transactions. At the moment of writing (August 2017), the largest blockchain they had analysed was Bitcoin. Bitcoin was at the time 478,559 blocks large representing 130GB of storage, but their analysis were scoped only to the 22GB of transactions. Ethereum contains, at the moment of writing this paper (January 30 2019), more than 7 million blocks representing 266GB of storage. It is not clear for us how Ethereum's constraints and its smart-contracts can easily fit in such a memory strategy. Bartoletti et al. [4] propose a general framework to access blockchain data previously stored in a database. In their approach, they export all blockchain data (around 500k blocks, 300GB of Bitcoin) into a NoSQL (mongoDB) and a relational (MySQL) database. Such imports take around 9 hours in each implementation, and querying them takes between 50 minutes and 3.5 hours, respectively. While they claim support for smart contracts and Ethereum, they did not report any benchmarks on it.

In prior work, we have also proposed the creation of a SQL-like query language to extract arbitrary data and analyse changes in user-defined smart contracts [START_REF] Bragagnolo | Ethereum Query Language[END_REF]. Like other researchers [START_REF] Hoang | Research Directions in Blockchain Data Management and Analytics[END_REF], we have identified how the performance problems of analysing the blockchain gets worse day after day. To perform fast searches by semantic attributes in such a query language, the usage of indexes is crucial. These indexes must be built by doing a full scan of all blocks, transactions, accounts and smart contracts in the blockchain, and building from them a smart data structure that can be quickly accessed. For instance, if we want to be able to access to all blocks that were created between two dates, an index data structure could be created to quickly discard all blocks that are outside the desired range. This kind of index data structure can be used for this and many other analyses, for example to access all the transactions that modified a given contract, or all the blocks where a given user/account received money. To build such indexes, however, a fast scanning infrastructure needs to be setup. Furthermore, such infrastructure needs to be kept up-to-date with the status of the blockchain. In our previous work we lacked such support.

BACKGROUND

Hardware advances in storage capacity and CPU processing have given rise to the concept of Big Data, characterized by the so-called 3 Vs (Volume, Velocity and Variety). As a result, novel software platforms have emerged to analyze and store such large data sets in a scalable way. The two most prominent programming models are Hadoop Map/Reduce [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF] and Apache Spark [START_REF]Apache Spark[END_REF], which typically embrace a batch-oriented data processing to achieve a high parallelisation of data analysis.

Current trends indicate that the volume, velocity and variety of data are increasing quickly due to an explosion on diversity and number of sources of information (as a result of the digitalization of data, e.g. smart objects and sensors, interconnectivity of data and popularity of social media data [START_REF] Cukier | Big Data: A Revolution That Will Transform How We Live, Work, and Think[END_REF]). In this paper we identify blockchain technology as a new type of sources of information in the Big Data revolution. Blockchain platforms like Ethereum contain a high amount of heterogeneous data growing at a high rate, presenting a clear example of the three Vs of BigData: Volume, Variety, Velocity. In this paper, we explore the usage of Big Data technology to enable blockchain analysis. A recent publication on research directions in blockchain analytics has also suggested the use of Big Data technology but, to the best of our knowledge, this is the first approach in realising such a vision [START_REF] Hoang | Research Directions in Blockchain Data Management and Analytics[END_REF]. In this section, we provide some details on the programming model we employed in this work, namely the Map/Reduce model [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF], and the programming platform in which we conducted this research.

The Map/Reduce Model

The Map/Reduce model can be used in a variety of different data analysis, and offers a simple API to developers. In this model, the developer structures her program in two functions: a map and a reduce. The map function transforms an initial set of values to an intermediate result, encoded in key/value pairs. The reduce function aggregates data to compute a final result. In particular, reduce transforms the intermediate results generated by map, grouped by their key, into to a final result for each of the keys. A Map/Reduce platform will take as input some collection of data, partition it and run the developer's map and reduce functions on each partition in different machines or processes.

As a concrete example, consider a grep which searches the lines containing a pattern in a file. The code snippet in Listing 1 implements such an operation on the Map/Reduce model. The map function checks if a particular line includes a pattern, and returns the line if the pattern is present, or NIL otherwise. The reduce function filters the mapped lines to remove the values that are NIL, returning all of the lines that contain the pattern. In this paper, we employed a Map/Reduce framework for Pharo Smalltalk called Port [START_REF] Marra | A debugging approach for Big Data applications in Pharo[END_REF]. Port defines a Smalltalk API to define and execute Map/Reduce programs deployed on top of the Apache Hadoop platform. In particular, developers need to write the map and reduce function in Smalltalk and then Port manages the execution of such a Map/Reduce application using master and worker nodes deployed on Hadoop Yarn [START_REF]Apache Hadoop YARN[END_REF], and implements all communication details between them. In this section we describe how a blockchain analysis application can be created using the Map/Reduce model. Figure 1 our Map/Reduce architecture which has been adapted to work with the Ethereum platform. The Map/Reduce framework provides a master node that splits the data and assigns map and reduce jobs to different worker nodes. We also setup a single database instance and a single Ethereum client, accessible by all these nodes in the same network. Each map worker queries the Ethereum platform to access blockchain data and sends it back to the master node. Then, each reduce worker receives many map results together, and performs a bulk insert in the SQL database. Once the data is stored in the database, we can use the information that was indexed (e.g., the hash of the block, or a particular property of the block or contract, ...) to query the blockchain and gather the requested information.

APPLYING BIG DATA TECHNOLOGY TO BLOCKCHAIN ANALYSIS

We will now illustrate the different parts of our Map/Reduce architecture to implement an indexation algorithm that scans the entire blockchain and stores a block index into the centralized SQL database.

An Indexation Algorithm

We have developed an indexation algorithm that uses a relational database to store indexed data. This choice was based in two reasons: to benefit from an existing and mature index implementation, and to focus on the problems of parallelism and high load. However, our overall architecture does not prevent alternative indexing data structures like a BST (Binary Search Trees) [START_REF] Bragagnolo | Ethereum Query Language[END_REF] to be implemented instead.

Our index has the structure of a relational table with standard database indexes. For example the table representing the block index has the block's hash, but also a timestamp and its parent block's hash as shown in Figure 2. The two latter columns are indexed, so we can do fast queries on blocks by both timestamp and their parent blocks.

To setup such an index, our core algorithm performs a full scan of the blockchain inserting all the corresponding values in our database. With a asynchronous I/O sequential scan approach, our algorithm takes about 30 milliseconds per indexed entry, meaning that the expected time for indexing the current load of Ethereum blocks, would take about 2 days and 11 hours.

To avoid a sequential scan, we decided to parallelize such processing using a Map/Reduce model as described below.

Parallelizing Blockchain Indexation

We now detail how we parallelized the indexation algorithm using a Map/Reduce model. This basically boils down to expressing the algorithm in terms of the map and reduce functions of the Map/Reduce model previously explained. Listing 2 illustrates our functions in pseudocode for the sake of clarity. In the next subsection, we detail the concrete implementation in Port. Such an approach allows us to parallelize the mapping of the indexes and to execute in parallel different database stores, minimizing the calls by executing one for each partition of data (and not one for each index).

Implementation

Figure 3 illustrates the main different components of our Map/Reduce architecture. The application is programmed in Pharo Smalltalk, on top on top of the Port Map/Reduce framework, which handles the management of master and worker nodes. Port provides us integration with the Map/Reduce environment offering an API to implement a Map/Reduce by the means of map and reduce methods in Smalltalk. Port uses Hadoop Yarn [START_REF]Apache Hadoop YARN[END_REF] to manage the allocation of master and worker nodes in a cluster.

We use Geth5 as blockchain data node. Internally, the communication with Geth is managed by the Fog Ethereum driver 6 .

Listing 3 illustrates our application implementation in Port. Our Map/Reduce application functions accept an additional parameter, an indexBuilder, that provides the database connection and all of the information about the properties that are being indexed. Such parameter is identical in all of the parallelized map and reduce functions. The map:parameter: (lines 1 to 6) obtains a block given a block index, obtains the indexed property from the given block, and returns a key value pair in the form (block index, indexed value). For example, in the case of indexing blocks by timestamp, the indexed value is the timestamp of the obtained block. The Differently from a classical Map/Reduce framework, Port does not enforce the initial data to be in a key-value store. This allows us to map directly on simple indexes, generating an intermediate key-value store where each key has only one value. It also does not always apply a group-by between map and reduce, as happens in classical Map/Reduce platforms. In our case this avoids an unneeded shuffling of data over the network, that would only slow down the blockchain analysis application.

EXPERIMENTS

In this paper, we aim to tackle the difficult process of indexing the full blockchain by applying a Map/Reduce algorithm on top of Port, a Map/Reduce framework for Pharo. We conducted different experiments to check whether our Map/Reduce architecture and implementation is a scalable architecture for blockchain analysis. A first experiment, called architecture, tries to find the best architecture (in the terms of amount of workers) to index the full blockchain. A second one, called analyzing the full blockchain, analyses how our approach scales when trying to analyze the full Ethereum blockchain.

Setup

We run our experiments on a cluster of one root node and ten identical slave nodes. Each slave node presents the following specifications:

-Processor: Intel Xeon CPU E3-1240 @ 3.50GHz (4 cores, 8 threads) -Ram: 32 GB -Storage: 200 GB SSD The root node has the same specification as the slave nodes, but it has enhanced storage. All the nodes are connected through a 1 GB/s local network.

In our configuration, the root node runs:

-A Postgres server that handles the database.

-An instance of Port to handle the external communication with the different nodes. One of the slave nodes runs exclusively Geth, the blockchain data node. The rest of the slave nodes ( 9) are used to deploy Pharo containers through Yarn. There is always one container running a Yarn application, one running a master, and a global maximum of 70 containers available to run a worker. The number of actual

Experiment 1: Architecture

Indexing the whole blockchain, even using a Map/Reduce approach, is not trivial. Not only does the performance of the indexing depends on the number of workers that can parallelize the work, but also on how these workers communicate with the blockchain data node (during the mapping phase) and with the database (during the reduce phase). While the employed database (i.e. Postgres) is designed for concurrent access and can scale well, the blockchain data node (i.e. Geth) allows us to query the blockchain using the RPC protocol, which can support only a limited amount of concurrent calls. Given our setup, Geth runs in a dedicated slave node with a total of 8 threads.

In this first experiment we test how the execution time changes when increasing the number of workers, while indexing a fixed amount of blocks (10.000 blocks). We increased the amount of workers from 4 to 48, 4 by 4. Each iteration is run 10 times, discarding the worst result.

Figure 4 shows the result of this experiment. The blue curve represents the average execution time of the full Map/Reduce, while the orange curve represents only the average time of the map phase. The black dashed line represents the averaged minimum amount of execution time (i.e. 43.6 seconds with 20 workers).

Looking at the total execution time, we observe that it quickly decreases from around 83 seconds (with 4 workers) to around 58 seconds for 8 workers. It then stabilizes between 47.7 seconds and 43.6 seconds when running with 12 to 48 workers, reaching its global minimum when running with 20 workers. As such, we conclude that 20 is the best amount of workers for this configuration.

Analyzing the relation between map time and total execution time, the graph shows that the two lines follow the same trend, being the map the biggest component of the whole execution time (between 83% and 92% of the total execution time). It is clear how the performance of the map is crucial to the performance of the whole indexation. We discuss more about the implications of this finding in Section 5.4.

Experiment 2: Analyzing the full blockchain

The results of experiment 1 show that the best amount of workers to use in our configuration is 20. In this second experiment, we use such configuration to check how our system scales to index the full blockchain. To do so, we gradually increase by a factor of 10 the number of blocks to index from 100 to the full Ethereum blockchain, that to the moment of our latest experiments (29/01/2019) consisted of 7.080.006 blocks.

Figure 5 shows how the execution time increases when increasing the number of analyzed blocks. The blue line represents the execution time, the grey line represents the execution time of the sequential implementation of the indexing algorithm. The values of the dashed gray line are proportionally projected, since the execution time would have been too high for experimenting it. The dashed black line represents what the proportional expected duration is, calculated on the result for 10.000 blocks. Both of the scales are logarithmic (loд 10 ), hence the two axis are proportional. The graph shows that the execution time grows proportionally, as the black dashed line, except for 100 and 1.000 blocks, where the execution time is heavily impacted by the overhead of Port.

Using 20 parallel single threaded workers, we managed to index the full Ethereum blockchain (to the moment the experiment was executed) in 7 hours 18 minutes and 47 seconds. While the experiments on up to 100.000 blocks were repeated at least 10 times, due to time limitations the experiments over 1 million blocks were executed only 4 times and the one on the full blockchain only once.

If we compare such results to a sequential approach (that we estimated would take more than 2 days), it is clear how the execution time of the indexation is lower using our Map/Reduce approach.

Discussion

The first experiment, Architecture aims to give some light on the limits of the Geth client and, therefore, in our first general architecture. The second experiment Analyzing the full blockchain, leverages the results of the first experiment for showing a good result on how a map reduce approach can improve the performance of blockchain analysis.

In the first experiment, we learned that the performance of the application does not scale to the number of workers, since between 12 and 48 workers the difference of execution time ranges in less than 10 seconds for 10.000 blocks. In fact, the map operation, which constitutes at least 83% of the execution time, performs calls to the blockchain data node (Geth), which does not handle optimally concurrent requests, since it is not designed to be used for analysis on the blockchain. This shows how making the analysis of the blockchain faster and scalable is not a mere question of parallelising, but also of finding the right configuration.

Since our first experiments we were sceptical about the responsiveness of Geth, the blockchain data node. To confirm quickly our hypothesis, we launched our algorithm with the maximum possible amount of workers deployed (70 worker nodes), one blockchain data node and one database node. Profiling the network performance exposed the limits of the blockchain data node towards the processing external requests, congestioning the network and overconsuming operating system resources. Hence, we conducted our first experiment starting from a little amount of nodes, to discover the limits of our architecture.

We believe that increasing the number of blockchain data nodes could further improve the map execution's time, heavily impacting the overall duration of the indexing. In our configuration, we found that the best ratio worker/blockchain data nodes is 20x, or, if we count on the physical threads the workers and the blockchain data nodes were using, around 2.5 workers for threads of the blockchain data node. However, further experiments are required to verify that such a ratio remains the same when increasing the number of blockchain data nodes (or the amount of cores it is using). In fact increasing such a number can increase the amount of network usage in the whole cluster, other than change the storage footprint of the blockchain that would need to be replicated in the different blockchain data nodes. We are planning to execute different experiments to investigate how such ratio can change, and what can be the impact on the overall execution time.

Overall, we believe that the experiments conducted on our approach shows how the use of a Big Data framework for parallelisation can improve the blockchain analysis.

CONCLUSION

In this paper, we explored the application of Big Data techniques to enable blockchain analytics. In particular, we applied Map/Reduce for indexing Ethereum blocks. To the best of our knowledge, this is the first work employing Big Data techniques to enable scalable blockchain analytics. Our experiments show that Map/Reduce can efficiently be used to index the blockchain, improving the run-time performances of the indexing process by 7.7 times with respect to a sequential implementation.

Our current architecture includes one node running a SQL database (i.e. Postgres), one node running the blockchain data node (i.e. Geth), one master node coordinating different worker nodes that run the map and reduce tasks. Our Map/Reduce application is written using Port, a Map/Reduce framework for Pharo Smalltalk. In such application, we parallelize the queries to the blockchain data node in the map function, trying to minimize the impact of such queries on the run-time. In the reduce function, we then perform a bulk insert of the indexed data in the database. By testing different setups of such architecture, we identified Geth as being a performance bottleneck, discovering the right ratio of worker nodes (that are consuming data from the blochckain) and blockchain data nodes.

Despite of the promising results shown in this paper, we believe it is too early to say that our architecture reached the limit of the Map/Reduce efficiency to enable blockchain analytics. Further research is needed to explore the best configuration for such an architecture. In particular, we are interested in exploring the following research directions:

• Exploring the advantages and limitations of using multiple data nodes. • Exploring the utility of Map/Reduce on the resolution of index based queries. • Exploring the limits of Map/Reduce as an online algorithm for continuous indexation.

Listing 1 :

 1 Pseudocode of a grep in Map/Reduce. 1 map ( fileLine , pattern ) { 2 if ( fileLine . contains ( pattern ) ) 3 then return fileLine 4 else return NIL } 5 reduce ( mappedLines ) { 6 return mappedLines . filter ( line => line ! = NIL ) }

Figure 1 :

 1 Figure 1: A Map/Reduce Architecture for Blockchain Analysis

Figure 3 :

 3 Figure 3: Our Map/Reduce Architecture

  reduce:parameter: (lines 8 to 9) receives a collection of pairs produced by the map:parameter: function and uses the builder to store those pairs in the corresponding database table.

Listing 3 :

 3 Pharo implementation of indexing algorithm 1 MRIndexingApp >> map : blockIndex parameter : aBuilder 2 MRIndexingApp >> reduce : pairs parameter : aBuilder 9 aBuilder storeIndexedValues : pairs

Figure 4 :

 4 Figure 4: Execution time of indexing 10.000 block increasing the number workers

Figure 5 :

 5 Figure 5: Execution time of indexing an increasing amount of blocks with 20 workers.

https://www.ethereum.org/

https://www.hyperledger.org/projects/fabric

https://nem.io/

https://wavesplatform.com/

https://github.com/ethereum/go-ethereum/wiki/geth

https://github.com/smartanvil/Fog

ACKNOWLEDGMENTS

Matteo Marra is funded by a SB PhD grant at FWO -Research Foundation -Flanders. Project number: 1S63418N.

We would like to thank UTOCAT for incentivating this work.


INTRODUCTION

The injection stretch-blow molding process is the most widely used in PET bottle production [START_REF] Rosato | Blow molding handbook[END_REF]. This process is segmented in two main stages. First, a cold preform is heated in an infrared oven above glass transition temperature (~80°C). Then, the softened preform is inflated and stretched with assistance of a rod in a mold which has the shape of the desired bottle.

A few heat transfer models in an infrared oven have been developed [START_REF] Lebaudy | Temperature distribution in poly(ethylene terephtalate) plate undergoing heat treatment, Diffusion influence : 1, Theoretical approach[END_REF][START_REF] Shelby | Effects of infrared lamp temperature and other variable on the reheat rate of PET[END_REF]. In 1993, DiRaddo et al [START_REF] Diraddo | Dynamic modeling of the preform reheating stage in injection blow molding[END_REF] proposed a 2D finite element model where the fraction of energy leaving the lamp surface and reaching at the preform surface (view factor) is estimated analytically. Hartwig [START_REF] Hartwig | Simulation of the stretch blow molding process and characterisation of t he relevant material behaviour[END_REF] included the influence of back and f ront reflectors and the displacement of the preform in the oven. Monteix carried out numerical simulations of the temperature distribution of preform using control-volum method and compared it with experimental data [START_REF] Monteix | Experimental study and numerical simulation of preform or sheet exposed to infrared radiative heating[END_REF].

The blowing step has been also the subject of significant researches within the last two decades [START_REF] Schmidt | Numerical simulation of polyester stretch-blow molding[END_REF][START_REF] Wang | 3D viscoplastic FEM simulation of a stretch blow molding process[END_REF]. But must of the papers use a hyperelastic behavior based on Mooney-Rivlin or neo-Hookean potentials. Queen's University of Belfast has programmed a special Abaqus ® version dedicated to the blow molding simulation using Mooney-Rivlin, Ogden or Buckley hyperelastic laws [START_REF] Menary | Modeling of Poly(ethylene terephtalate) in injection stretch blow molding[END_REF]. Recently, Gorlier [START_REF] Gorlier | Experimental and theoretical study of uniaxial deformation of amorphous P.E.T. above glass transition temperature[END_REF] or Marco [START_REF] Marco | Caractérisation multi-axiale du comportement et de la micro -structure d'un semi-cristallin : application au cas d u P[END_REF] pointed out that the Edwards Vilgis potential [START_REF] Edwards | The effect of entanglements in rubber elasticity[END_REF] allows to capture the deformation of PET subjected to multi-axial loads. Moreover, the thin shell models reach their limits regarding the more and more complex bottle shapes.

As a consequence, the two steps of the stretch blow molding process will be computed in this work with a 3D finite element model implemented in the Forge3 ® commercial software. First, a modeling approach capable of predicting the threedimensional transient temperature profile in the preform during the conditioning step is presented, taking into account conduction, convection and radiation. The mechanical problem is considered using Newtonian or hyperelastic Mooney-Rivlin laws implemented with a velocity/pressure formulation.

INFRARED HEATING MODEL

In its classical form, the heat balance equation can be written :

( ) ( ) r p q T k dt dT c ⋅ ∇ - ∇ ⋅ ∇ = ρ (1)
Notations are detailed in the nomenclature at the end of the paper. To calculate the density radiative heat flux, a ray tracing method coupled with a Beer-Lambert approximation has been implemented. It allows to calculate with accuracy the absorption of the radiation inside the material. Along each ray emitted, the Beer-Lambert law can be applied :

d e q q l A r - = 0 (2)
ABSTRACT: In this work, both the heating stage and the blowing stage of the blow molding process are numerically modeled. The heat transfer between the infrared oven and the preform is modeled using a ray tracing method. The cooling fan effect is taken into account thanks to a forced convection coefficient.

Regarding the blowing step, a Mooney-Rivlin hyperelastic model has been implemented in Forge3 ® software in order to account for the rheological behavior of the polyester. The finite element model is based on velocity pressure formulation and tetrahedral elements. In order to validate the implementation of the hyperelastic behavior, computations are compared to Mooney-Rivlin analytical model results for tube free inflation. Finally, the global blow molding process of a PET bottle is studied.

Key words: Blow molding modeling, hyperelasticity, radiation, ray tracing.

3D finite element modeling of the blow molding process

C. Champin 1,2 , M. Bellet 1 , F.M. Schmidt 2 , J.-F. Agassant 1 , Y. Le Maoult 2
The intensity of radiation emitted by a lamp at 2415K is about L = 1.66x10 5 Wsr -1 m -2 [START_REF]Training period report[END_REF]. As a consequence, the emissive power of radiation for each beam can be calculated :

Elt W beamsS Nbr LS q π = 0 (3)
The energy absorbed inside the material depends on the average absorption coefficient which can be determined using the Planck function and the measured value of the spectral absorption coefficient ( )

( ) ∫ ∫ ∆ ∆ = λ λ λ λ λ λ λ d T L d T L A A 0 0 (4)
The halogen lamp is modeled by its tungsten filament and its ceramic reflector. For each surface node, a constant number of emitted rays is defined. The orientation and the energy of each ray is defined in order to get an isotropic radiation in the half space determined by the outward normal vector.

In a first step, to calculate the interception between a ray and the preform heated, a three equations with three unknowns system is solved considering the beam equation and the boundary element area of the preform. If the beam intercepts f irst the back reflector, its direction is recalculated considering a specular reflexion and tested again with the object. In a second step, the ray's trajectory and the heat flux absorbed are defined using the connectivity of the mesh and the Beer-Lambert approximation.

The irradiance calculation has been validated comparing the Forge3 ® results with an analytical model of an element on plane facing a cylinder [START_REF] Leuenberger | Compilation of radiation shape factors for cylindrical assemblies[END_REF]. As shown in figure 1, the agreement is fair (about 10%). The calculation of the amount of energy received by the plate is validated using Philips measurement [START_REF]Training period report[END_REF]. Like the comparison with the analytical model, the final error is about 10% (figure 2). Finally, the temperature evolution in the preform heated with one lamp (figure 3) or with six lamps (figure 4) is obtained, taking into account the natural convection and the ceramic reflector of the lamp.

MECHANICAL MODEL

As shown by many authors [START_REF] Gorlier | Experimental and theoretical study of uniaxial deformation of amorphous P.E.T. above glass transition temperature[END_REF][START_REF] Marco | Caractérisation multi-axiale du comportement et de la micro -structure d'un semi-cristallin : application au cas d u P[END_REF], the PET exhibits a hyperelastic rheological behavior above its glass transition temperature. Assuming incompressibility and Mooney-Rivlin potential, the expression of the Cauchy stress tensor is [START_REF] Germain | Cours de Mécanique des Milieux Continus[END_REF] :

1 ' ' ' 2 2 - - + - = + - = B C B C I p I p α σ σ ( 5
)
This rheological model has been implemented in Forge3 ® using a velocity/pressure updated Lagrangian formulation together with P1+/P1 tetrahedral finite elements.

This results in a set of non linear equations, the unknown of which are the nodal velocity field and the nodal pressure field. This set of equations is solved by a Newton-Raphson method. For validation purpose, an analytical solution has been developed [START_REF] Champin | [END_REF] in an axi-symmetric thick tube inflated with an imposed internal pressure. The lower and upper sections of the tube are supposed to slip freely along fixed horizontal planes. The polymeric tube is assumed to have a Mooney-Rivlin behavior. In figure 5, both numerical and analytical internal radius versus pressure are plotted as well as error.

The agreement is fair (maximum error 4%).

BLOWING MODEL

An isothermal free inflation simulation of a preform has been performed using a neo-Hookean behavior ( 0 = α in equation ( 5)). The parameters used for this calculation have been referenced in table 1. Moreover, the relaxation phase due to the hyperelastic behavior can be checked numerically : if the pressure applied inside the preform decreases to zero, then, the preform relaxes to its initial shape.

As a second example, a confined blow molding simulation (without rod) has been done. We assume an isothermal Newtonian law behavior. The mould used is a prototype one developed at CROMeP (table 2). It produces 50 cl water bottle. As for the tube inflation, just one eighth of the preform has been meshed (figure 7).

CONCLUSION AND PROSPECTS

A numerical thermo-mechanical model has been developed using Forge3 ® velocity/pressure formulation and tetrahedral elements. From a thermal point of view, a ray tracing method coupled with a Beer-Lambert approximation has been implemented. This method has been compared with an analytical model of view factor calculations and irradiance measurement on a plate. Ceramic reflector and tungsten filament have been modeled. The calculations confirm the major role of the back reflector on each lamp. From a mechanical point of view, a Mooney-Rivlin hyperelastic behavior law has been implemented with a finite element formulation including consistent tangent stiffness terms. The model has been validated by comparison with an analytical calculation of tube inflation. Future work will consist in testing alternative hyperelastic potentials and modeling the stretch step in order to compare the final thickness distribution with experimental data.

Fig. 1 .

 1 Fig. 1. Irradiance of a plate (350 x 350 mm), 100mm distance from a lamp (350 mm length / 0.81 mm diameter).

Fig. 2 .

 2 Fig. 2. Irradiance of a plate (600 x 600 mm), 100 mm distance from a lamp (280 mm length / 2.15 mm diameter)

Fig. 3 .

 3 Fig. 3. Surface temperature versus time at mid -height of the preform depending on number of beams and back reflector.

Fig. 4 .

 4 Fig. 4. Surface temperature versus time at mid -height of the preform depending on number of lamps.

Fig. 5 .

 5 Fig. 5. Internal radius variation versus pressure.

Table1.

  Figure 6 presents successive steps of the free blow molding process.

Fig. 6 .

 6 Fig.6. Successive configurations during free inflation.

  Length of the print = 119 mm

Fig. 7 .

 7 Fig. 7. Geometries of the preform and the mold.

Figure 8

 8 Figure8shows the inflation of the preform as a function of time and especially the development of the contact between polymer and mold.

Fig. 8 .

 8 Fig. 8. Successive configurations during blow molding.

kqx

  Thermal conductivity (W m -1 K -1 ) r q Density radiative heat flux (W m -2 ) 0 Initial density radiative heat flux value (W m -2 ) A Average absorption coefficient (m -1 ) l Distance covered inside the material (m) Spatial coordinate (m) i v i-th node velocity (m s -1 ) P ∆ Pressure applied (Pa) R Internal radius (m) S External radius (m)

Ecole des Mines de Paris -CEMEF -BP207 F-06904 SOPHIA ANTIPOLIS Cedex e-mail: cedric.champin@ensmp.fr; michel.bellet@ensmp.fr; jean-francois.agassant@ensmp.fr

Ecole des Mines d'Albi -CROMeP -F-81013 ALBI Cedex 09 e-mail: fabrice.schmidt@enstimac.fr; yannick.lemaoult@enstimac.fr

ACKNOWLEDGEMENT

This study has been supported by Nestlé Waters Management & Technology (PTC Vittel). Special thanks go to G. Denis, K. Hartwig and C. Boulay for their help in this work.


Maxime Bordival

Sylvia Andrieu

Fabrice M Schmidt
email: schmidt@enstimac.fr

Yannick Le Maoult
email: lemaoult@enstimac.frbordival@enstimac.fr

Serge Monteix
email: monteix@philips.com

Optimization of infrared heating system for the thermoforming process

Keywords: thermoforming process, infrared heating, control volume method, raytracing method, Rosseland model, SQP optimization thermoforming process, infrared lamp, radiative heat transfer, control-volume method, raytracing method, Rosseland model, SQP optimisation

INTRODUCTION

In a typical thermoforming process, a polymeric sheet is heated to its forming temperature and brought into contact with a mold of the desired shape. Owing to the low polymer thermal conductivity (0.29 W/m/K for PET), the heating step is performed using infrared lamps, which enable short thermoplastic sheet heating times. The final thickness of the part, and consequently the final product quality, is drastically controlled by the initial temperature distribution inside the sheet. For this reason, optimizing the heating stage is crucial in order to ensure an improved process control. The thermal problem during thermoforming process results from a coupling of conductive, convective and radiative three-dimensional unsteady transfers. Moreover, the modelling of heating device, composed by different rows of lamps associated to reflectors, is complex. Thus, to find optimum process parameters, inverse methods are generally used.

Thermoforming optimization approaches

Different recent optimization process investigations have been studied by using commercial software such as POLYFLOW ® [START_REF] Monteix | Optimization of the Thermoforming Process: A Few Industrial Examples[END_REF] and T-SIM ® [START_REF] Novotny | Optimization of Thermoforming[END_REF]. Some investigators, as Wang and Nield [START_REF] Wang | Temperature Optimization for improved Thickness Control in Thermoforming[END_REF], developed their own thermoforming simulation software. Both of these studies, based on mechanical stage simulation, focused on evaluating the possibility to obtain uniform thickness by controlling the initial polymer sheet temperature field. However, ovens parameters were not linked to the temperature distribution. Duarte and Covas [START_REF] Duarte | Heating Thermoplastic sheets for Thermoforming Solution to Direct and Inverse Problems[END_REF], and Throne [START_REF] Throne | Radiant Heat Transfer in Thermoforming[END_REF], proposed algorithms to solve the inverse heating problem, in order to find the heater temperature pattern associated to a uniform sheet temperature field. Both of these studies were based on solving the onedimensional unsteady heat conduction problem, using radiative and convective boundary conditions. As a conclusion, no inverse method has been used to find oven optimum geometric parameters.

ABSTRACT: Determining an optimum processing window in thermoforming process is critical in the aim of achieving high quality parts. The infrared heating step is crucial, the final thickness distribution of the thermoformed part being closely related to the initial temperature field inside the sheet. To reduce the time needed to dimension the heating device, an automatic optimization method of ovens geometric parameters has been developed. In a first time, a simple analytical model, coupled to a non linear constrain optimization method (Sequential Quadratic Programming), permits to find the best set of parameters, according to a cost function representing for example the heat flux uniformity. Then, with these optimized parameters, an accurate raytracing method is used to compute the irradiation resulting from the interaction between lamps and thermoplastic sheet. Finally, control volume method is implemented to solve the three dimensional transient heat transfer equation, where the radiation source is approximated by a diffusion Rosseland model. ABSTRACT: Determining an optimum processing window in thermoforming process is critical in the aim of achieving high quality parts. The infrared heating step is crucial, the final thickness distribution of the thermoformed part being closely related to the initial temperature field inside the sheet. To reduce the time needed to design the heating device, an automatic optimization method of ovens geometric parameters has been developed. In a first time, an analytical model, coupled to a non linear constrain optimization method (Sequential Quadratic Programming), permits to find the best set of parameters, according to a cost function representing for example the heat flux uniformity. Then, with these optimized parameters, an accurate raytracing method is used to compute the irradiation resulting from the interaction between lamps and thermoplastic sheet. Finally, control-volume method is implemented to solve the three-dimensional transient heat transfer equation, where the radiation absorption is approximated by the diffusion Rosseland model.

Goals of the present study

We propose an inverse method, permitting to optimize the oven geometry, and to obtain a threedimensional temperature pattern at the end of the heating stage. This method can be decomposed into three coupling steps:

• Optimization of the heat flux over the sheet surface, using an analytical model coupled to a non linear constrain optimization method (Sequential Quadratic Programming). The goal is to find the best set of oven geometric parameters. • Accurate simulation of radiative transfers, based on raytracing method, with optimized parameters obtained in the first step. • Computation of the 3D temperature pattern, using THERMORAY, a heat transfers solver (developed in CROMeP [START_REF] Monteix | Experimental Study and Numerical Simulation of Preform or Sheet Exposed to Infrared Radiative Heating[END_REF]) based on control-volume method. The absorption of the irradiation calculated in the second step, is approximated using the Rosseland model.

HEAT BALANCE EQUATION

In its classical form, the heat balance equation can be written as:

) .( ) .( ) .( ) .( → → → → ∇ - ∇ ∇ = ∇ - ∇ = r r c p q T k q q dt dT c ρ (1)
where T = temperature, q c = conduction heat flux, q r = radiation heat flux, ρ = density, c p = specific heat, k = thermal conductivity. To compute the temperature distribution through the thickness of the sheet, this equation is solved using a 3D control-volume software, named THERMORAY [START_REF] Monteix | Experimental Study and Numerical Simulation of Preform or Sheet Exposed to Infrared Radiative Heating[END_REF].

Control-volume method

The sheet is meshed using cubic or hexahedral elements called control volumes. The equation ( 1) is integrated over each control volume and over the time from t to t+∆t:

dt d n q dt d n q dt d t T c t e r t e c t e p Γ - Γ + = Ω ∂ ∂ ∆ Γ → → ∆ Γ → → ∆ Ω ) . ( ) . ( ρ (2)
where Ω e = control volume. The unknown temperatures are computed at the cell centres of each element. Different ways exist to approximate the radiation heat flux absorption. One of them is the Rosseland model.

Rosseland approximation

This method allows to approximate the irradiation absorption. Though it avoids long computation times, it is reserved to optically thick media [START_REF] Modest | Radiative Heat Transfer[END_REF]. The main idea of this approximation is to consider that the irradiation is diffused inside the medium. Thus, the heat balance equation is solved in a pure conduction approach; in other words, the radiation heat flux reduces to:

∇ ∇ - = ∇ → → T T k q Ross r ) ( . ) ( . (3)
Optical properties inside the polymer are taken into account via the Rosseland conductivity defined as:

3 ) ( 3 16 T T K k Ross Ross σ = (4)
where σ = Stefan-Boltzmann constant, K ross = Rosseland mean coefficient obtained versus temperature by integration over the frequency:

ν ν ν ν ν κ d dT T dB d dT T dB T K Ross ∞ ∞ × = 0 0 ) , ( 1 ) , ( ) ( 1 ) ( 1 (5)
where κ = absorption coefficient, ν = frequency, B = black body intensity (Planck function).

Raytracing method

In order to compute the irradiance, a raytracing method has been implemented [START_REF] Cai | Development and application of Radiative Heat Transfer using raytracing method[END_REF]. This method permits to model accurately the IR lamps, taking into account emitters geometry and optical properties, which depend on temperature and frequency, as well as back reflectors.

Analytical model

Raytracing method is very accurate. Unfortunately, it is also very slow. In order to reduce the time consumed by the optimisation step, an analytical model is used. It is based on the expression of the irradiance of a surface element, from an isotropic emitter point. Meshing Nl lamps in Ns elements, and using irradiance additive property, the irradiance E on a rectangular sheet element is given by [START_REF] Monteix | Experimental Study and Numerical Simulation of Preform or Sheet Exposed to Infrared Radiative Heating[END_REF].

[ ]

+ = L S N i N j ij ij ij ij x H H P E 2 3 2 2 4 1 π (6)
Where P ij = power emitted by each source point, and H ij , x ij defined as shown in figure 1.

Fig. 1. Analytical model: adopted notations.

INFRARED OVEN OPTIMIZATION

The analytical computation of irradiance is coupled to an optimization method, to automatically modify oven geometric parameters at each iteration, as shown in figure 2. The problem being non linear, continuous, and constrained, an adapted Sequential Quadratic Programming method is used in Matlab ® [START_REF]Optimization Toolbox For Use With Matlab®[END_REF].

Cost function

The optimisation objective is mathematically represented by a cost function F to minimize, function of oven geometric parameters X. We propose to establish the function ( 7) related to radiation heat flux uniformity at the surface of the sheet, and to a desired mean value of irradiance.

OBJ OBJ E E X E X X F ) ) ( ( 5 . 0 ) ( 5 . 0 ) ( 0 - + = σ σ (7) 2 , ) ( 1 )
( - = n j i ij E E n X with σ (8)
where E(X) = mean irradiance computed at each iteration, E OBJ = desired mean irradiance, σ(X) = standard deviation of computed irradiances at each iteration, σ 0 = standard deviation of computed irradiances with initial conditions, n = number of sheet surface elements. This function permits to obtain not only a uniform irradiance field, but a desired mean value of irradiance too.

Parameters

Both process and geometric parameters can be included in the optimisation. In this study, we focus on geometry only. Parameters concerned are: IR lamps positions for the two rows (X s , Y s , X i , Y i ), the distance between lamps and sheet (H), and the distance between lamps (D), as shown in figure 3.

Constraints

Geometric parameters are bounded. For example the distance between lamps must be greater than 20 mm, not to damage IR lamps. These constraints can be expressed as inequality relations between several parameters. For example: X s + 3D < 550.

APPLICATION

The method developed is used to optimize the geometry of an oven constituted of eight lamps Philips 1 KW clear, according to the cost function [START_REF] Modest | Radiative Heat Transfer[END_REF]. For the raytracing, 563 million of rays are followed. The value of the convective coefficient, applied at sheet surface during heating stage, is 15 W/m²K. The value of the PSB thermal diffusivity is 1.92 .10 -7 m²/s. As it is shown in figures 4, a more homogeneous temperature distribution over the sheet surface is obtained using updated oven geometry, even after 15 s heating. The maximum temperature difference is around 3°C.

CONCLUSIONS

A complete numerical model of the heating stage has been developed. Using THERMORAY, a control-volume software, and the raytracing method, the 3D heat balance equation is solved with convective and radiative boundary conditions. The irradiation absorption is approximated by the Rosseland model. Concerning the optimization, the SQP algorithm coupled to an analytical model, permits to design optimally ovens geometry, in order to increase the surface temperature uniformity.

The future work will consist in taking into account reflectors in optimization, which are generally used to increase the process efficacy. Moreover it is crucial to enlarge the model to optically thin media, and thus, to compute the irradiation absorption using a more accurate method.

Fig. 2 .

 2 Fig. 2. Heat transfer simulation / optimization coupling.

Fig. 3 .

 3 Fig. 3. Adopted configuration / parameters notation.

Fig. 4 .

 4 Fig. 4. Front surface temperature distribution after 15 s heating. Initial and updated oven geometry.

Fig. 5 .

 5 Fig. 5. Convergence curve.

  The objective mean value of irradiance is 6000 W/m². Dimensions of the PSB (PolyStyrene Black) sheet are: 350x350x2 mm. It is meshed into 1225 rectangular elements. Results are reported in table 1.

	Table1. Initial and final parameters
	Parameters (mm)	H	D	Xs	Ys	Xi	Yi
	Initial	100	100	-90	300	-25	-115
	Final	146	151	-88	Isotropic emitter 477 -23	-121
		d		Hij
			Xij
		dS	Non participative medium
	CPU time: 385 s (AMD Athlon XP 2600+, RAM: 500 Mo).
	Number of iterations: 26.
		Initial Process Parameters		Initial Geometry		(Xs,Ys)	D	Polymeric
								sheet
			Sheet and Infrared Lamps Meshing
		Analytical Model					IR lamp
	Irradiance field over the Sheet
					Updated Geometry
		Cost Function Computation				(Xi,Yi)
				Updated Process
				Parameters
	Sequentiel Quadratic
		Programming
		Optimum ?	No
		Yes			Raytracing / control-volume method
		Optimized
		parameters
					Temperature field

CROMeP -Ecole des Mines d'Albi CarmauxCampus Jarlard -81000 Albi, France URL: www.enstimac.fr e-mail: bordival@enstimac.fr; andrieu@enstimac.fr schmidt@enstimac.fr; lemaoult@enstimac.fr

Philips Lighting Chem. Montrichard -54700 Pont à Mousson, France URL: www.infrared.philips.com e-mail: serge.monteix@philips.com

ACKNOWLEDGEMENTS

This study was funded and supported by Philips, Lighting Division.


N Lebaal

S Puissant

F M Schmidt

Rheological parameters identification using In-situ experimental data of a flat die extrusion

Keywords: MLS, Numerical Simulation, Identification of rheologicals parameters, Extrusion

.

Global relative error between these both pressure is the objective function. This objective function is minimised by varying the rheological parameters. For this minimization two methods are used i.e. local response surface and global response surface. The rheological parameters obtained by these methods allow calculating the viscosity, validated by comparing an experimental viscosity measured on a capillary rheometer.

Introduction

The defects of extrusion (like the weld-lines, the bad distribution, a fairly uniform exit velocity distribution throughout the extrusion, problems of stagnation zones) are influenced by the geometry of the die of extrusion, the operating conditions such as temperature of regulation, flow rate and the rheological parameters of melt. In order to eliminate these problems we propose to determine a sufficiently robust method of optimization to be applied in industry. In a first step, in order to validate the method of optimization, we will identify the rheological parameters (behavior of the pseudo plastic melt) of the plastic melts starting from the points of measurement obtained in flat die.

A standard extruder at various flow rates feeds the flat die and two series of tests were carried out, with two thicknesses of die to cover a sufficient range to make rheological characterizations.

In spite of various utilities of plastics, the knowledge and the control of their behaviours remains always approximate.

Physicists tried to characterize the properties of these materials, thus wanting to reproduce their behaviours using a simple or complex model. A number of viscosity Corresponding author. E-mail address : lebaal@enstimac.fr models have been published, such as the well-known Power-law model, Ostwald-Waele model, Cross model, and Carreau model, and so on. These models are established on the basis of some types of fluid. Ji-Zhao Liang [START_REF] Liang | Characteristics of melt shear viscosity during extrusion of polymers[END_REF] he investigated the melt viscosity in steady shear flow of several polymers by employing a capillary rheometer and the characteristics of shear viscosity. K. Geiger and H. Kuhnle [START_REF] Geiger | Analytische Berechnung einfacher Scherstromungen aufgrund eines Fliebgesetzes vom Carreauschen Typ[END_REF] used a empirical correlation between density, pressure and temperature to identify a rheological parameter.

Our objective is to identify the rheology of a plastic melt directly from on-line production because that saves us much time and money, because in this way we will not be obliged to make measurement on dies capillary standard.

For that we will identify the rheological parameters, using algorithm of optimization named as response surface method, with the moving least squares approximate [5 -6], and can obtain approximation by using the design of experiments [START_REF] Goupy | Plans d'expériences pour surfaces de réponse[END_REF].

The optimization benchmark

The flat die is equipped with 4 pressure transducers, spaced 100 mm apart; however only the pressure difference between the first and last transducer was considered for the analysis. A first series of tests has been performed with the flat die. Measurements were made with two different slit heights at different extrusion speeds for one type of material (LDPE). The experimental data is used as listed in table 1 Fig. 1. Geometry of the Flat die the geometry of the die is W = 100 mm, L = 300mm. L total =500 mm. The thickness (h) of the die remains constant for a test, but it can be changed by taking in to account the correction factor (h/L) [START_REF] Agassant | Polymer rocessing. Principles and Modeling[END_REF].

Numerical simulation parameters

The measurement of viscosity cannot be obtained directly because the shear rate is unknown and the three parameters for the power-law. On the other hand, three physical parameters related to viscosity such as pressure, temperature and the flow rate are directly measurable. We will use a one-dimensional calculation programme by sections of finite difference method [START_REF] Agassant | Polymer rocessing. Principles and Modeling[END_REF], to calculate the variation of pressure with power-law for each point of measurement. Choice of this geometry by taking in to account the Correction factor[1] 0.1 and 0.05 as a function of the thickness to width ratio, h/W.

Design variables and objective functions

The power-law is defined as:

1 m T 1 0 e K (1)
There are three variables (K, m and ß) for a power law with a thermal dependence of the Arrhenius type.

The objective function is the difference between the measured pressures and calculated pressures obtained by simulated program of finite difference method.

The objective of this simulation is to determine the rheological parameters K, m and ß which makes it possible to minimize the difference between the calculated pressures and the experimental pressures. We chose as variable optimization the rheology of the polymer and the objective function represents the value of the sum of the variations of the calculated pressures and experiment pressures for each point of measurement.

J=(

2 exp npm cal exp ) p / ) p p ((
).

(2)

Optimization procedure

Choice of algorithm

The algorithm of optimisation must be carefully chosen when one single analysis requires several hours of CPU time. Nondeterministic or stochastic methods such as Monte-carlo method and genetic algorithm [START_REF] Renders | Algorithmes génétiques et réseaux de neurones[END_REF] can obtain global minimum but they need a lot of evaluations for the functions to converge. Decent methods require the computations, the gradients of the functions, such as BFGS [START_REF] Morales | A Numerical study of limited memory BFGS methods[END_REF] with ought constraints, and SQP [START_REF] Horowitz | Quadratic programming solver for structural optimisation using SQP algorithm[END_REF] based on Lagrangian methods, it can be shown that the solution of the nonlinear equality constrained optimization problem. The computation of gradients by finite difference is time consuming and depends on the perturbed parameters. For the above reasons we decided to consider a response surface method.

Response surface method

The method of response surface consists the construction of an approximate expression of objective function starting from a limited number of evaluations of the real function.

The main idea is to approximate the objective function through a response surface. To obtain a good approximation we used a moving least square method. In this method, the approximation is computed by using the evaluation points by design of experiments around the locus where the value of the function is needed.

Moving least-squares approximation

We will use moving least-squart interpolations that are used to approximate response surface of objective function. Let the dependent variable be ) x ( f in the domain S and the approximation of

) x ( f ~. P1 h W P2 L L total ) x ( a ) x ( p ) x ( f ~T (3)
The nodes are defined by x 1 ,…..xn where x I =(x I ,y I ) in 2D, x I =(x I ,y I ,z I ) in 3D; Where ) x ( p is a polynomial basis and ) x ( a T is the vector of coefficients. The polynomial basis of order 2 in two dimensions is given by

2 2 1 2 2 1 i ) y ( ) y )( x ( ) x ( ) y ( ) x ( 1 P (4)
In the moving least-square interpolation, at each point x, a is chosen to minimize the weighted residual:

2 T i i )) x ( a ) x ( P f ( w ) a ( J ) Pa F ( ) Pa F ( ) a ( J 2 1 (5)
wher i w is a weight function, such that w is non-zero over a limited support or domain of influence and is nonincreasing over. To find the coefficients, we obtain the extremum of J(a) by

0 ) a ( J ' So we have Z . B Q ) x ( a 1 (6)
where WP P Q T [START_REF] Renders | Algorithmes génétiques et réseaux de neurones[END_REF] and WQ P B T [START_REF] Horowitz | Quadratic programming solver for structural optimisation using SQP algorithm[END_REF] Z is the dependent variable.

Local response surface:

This method is based on traditional minimization by the algorithm of fig2, with use of a local response surface of variable in each point. We fixed a step of grid at 0.01* X 0 .

X 0 is the initial value of the variable of optimization The initial values are: K 0 =12000, m 0 =0.641370 and ß 0 =1078 °K, we used a uniform local grid i.e. "Design of experiments"

The constraints for the variables of optimization K, m and are respectively, 400 < K < 20000, 0.01 < m <1 and 500 < ß < 3000 °K for the first constraint, and n*P for the second constraint, with P is the step of the grid and n is the number of points of grid.

Fig. 2. Flowchart of the optimisation

The iterative process stops when the successive points are superpose with a certain tolerance.

Global response surface

The objective of this technique is to search for the global minimum on the field domain, the constraints is similar to the first constraints in local response surface method. Initially, in first iteration we constructed a grid of 5*5*5 points (125 points) for a global approximation, then to each iteration we constructed grid of 3*3*3 (27 pints) for the local approximation.

Results

Result of local response surface:

Figure 3 represents the evolution of the objective function and figure (4), ( 5) and ( 6) represents respectively the evolution of the coefficients of the law of rheology K, m and during the process of optimization. The optimal solution, which represents the minimal variation between the pressures obtained by finite difference calculation and the measured pressures,

In iterations 20, the value of the function objective obtained is 8.e-5 , with K= 9826.14, m= 0.4179 and ß = 1553.41, a time CPU is 1hours 42minutes, on a machine Pentium IV, 2.4 GHz, 512 Mo RAM, the results of the pressures are presented in table 2.

Optimal solution

Yes

Result of global response surface:

The solution is obtained after 6 iteration with a precision of 8.4e-3 and a time computing CPU of 1hours 13minutes on a machine Pentium IV, 2.4 GHz, 512 Mo RAM. We remarked in the table (2) that the solution obtained by the method (1) local response surfaces is more precise but with a more increased computing time where the difference in time between both methods is 30 min. On the other hand the global response surface method (2) avoided local minimum and we can also improve by varying the design of experiments and the method of approximation to optimize the computing time CPU. The figures ( 7), (13) illustrates the difference between the measured pressures and calculated pressures obtained by rheogical parameters. Which are identified by optimization using response surface method. It is noticed that there is almost no variation, which means that our method is very robust and gives very good results for the two approximations global and local.

Viscosity

We remarked well in the graph of viscosity obtained by the coefficients found by response surface method fig ( 8), ( 14), that it does not have an error between the points of measurements obtained in capillary rheometer and viscosity obtained by the method of optimization except in the zone that has high shear rate. As we have no knowledge of the inlet temperature so we used a simplification where we gave the same temperature of exit as of inlet temperature, on the other hand at this level of shear rate the polymer is reheated and this simplification is not completely right, but the variation between viscosities is always less important.

Concluding remarks

We can have the conclusion that since we obtained good results by the response surfaces method (global and local), so this allows us to know the rheology of a plastic directly in series of production and that saves us much time and money, because this way we will not be obliged to make measurement on standard capillary dies. s -1

Pa . s Pressures bar

Pressures bar

With this method we can find out the rheology of polymer starting from the true geometry even if this geometry is complex, by measurements of pressures, flows rate and the temperatures, we can determine the rheology of the material used.

Prospects

We have tried another method of approximation for approximating well the global response surface using the Krigeage method, which made it possible to approximate the global objective function, and saves us much computing time.

Using response surface method to optimize the geometry of die extrusion using calculate program 3D " REM3D ® "

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 7 .Fig. 8 .

 34578 Fig.3. Convergence of the objective function (result for the local method)

  Fig 9 represents the evolution of the function objective and fig 10, 11 and 12 represent respectively the variables of optimization K, m and . The optimum parameters are: K=10213, m= 0.4064, ß=1386.2617 °K.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Convergence of the objective function (results for the global method)

  .

	Table 1.
	Experimental data slit 5 and 10mm
	Q [Kg/h]	T [°C]	? P[bar] ? [Kg/m 3 ] h [mm]
	11,2	166	53	779	5
	60	175	99	779	5
	143	192	125	773	5
	239	210	138	767	5
	10,9	185	12,7	779	10
	69,2	195	26	772	10
	189,5	210	34,3	761	10
	339,5	229	36,7	755	10
	Conductivity k=0.115 Watt/°C m
	Specific heat Cp=1.90e03 K joules/°C m 3

Table 2 .

 2

	Results of optimization
	h [mm]	Q [Kg/m 3 ]	T [°C]	? P [bar]	? P[bar] local (1)	Residue n (1)	? P[bar] global (2)	Residue n (2)
	5	3,9e-6 166 53 52.11 2.8e-4 49.32 4.8e -3
	5	2,1e-5 175 99 96.69 5.4e-4 91.03 6.4e -3
	5	5,1e-5 192 125 121.6 7.3e-4 116.8 4.3e -3
	5	8,6e-5 210 138 132.9 1.3e-3 130.4	3e -3
	10 3,8e-6 185 12.7 12.49 2.7e-4 12.48	3e -4
	10 2,4e-5 195 26 25.15 1.e-3	24.89 1.8 e -3
	10 6,9e-5 210 34 34.69 1.2e-4 34.66 3.7 e -4
	10 1,2e-4 229 36.7 39.27 4.9e-3	40	8.1 e-3
	10.1. Pressures

Acknowledgments

This document was created with Win2PDF available at http://www.daneprairie.com. The unregistered version of Win2PDF is for evaluation or non-commercial use only.


G Dusserre
email: gilles.dusserre@enstimac.fr

F Schmidt

G Dour

G Bernhart

Thermo-mechanical stresses in cast steels dies during glass pressing process

Keywords: Glass Pressing Process, Numerical Simulation, Thermo-mechanical Stresses, Rosseland Diffusion Approximation

Cast steels dies used in glass pressing process present dispersal in their lifetime in spite of apparently similar pressing parameters. These dispersal does not seems to be due to casting defects because every mold presents such defects. In order to know if undetectable variations of any process parameter could produce such lifetime dispersal, a finite element simulation of glass pressing process is here proposed. In this paper, glass is supposed to be a Newtonian fluid, cast steel is considered as a thermo-elastic material, thermal transfer between glass and mold depend only on time and radiative heat transfer in glass is taken into account by Rosseland diffusion approximation. The validation of the results of this simulation is made by comparison of the calculated pressing force to the experimental one. Although the experimental data introduced in the simulation does not have been measured in the industrial process, the order of magnitude of the calculated pressing force is similar to the real one.

Introduction 

Glass pressing is a forming process that realizes hollow parts in only one stage (Fig. 1). The glass is heated up to 1500°C in a furnace and then cooled near 1250°C in the feeder. The gob is formed and fall into the lower mould preheated at 450°C. The gob shape then between the lower and the upper mould, activated by an hydraulic press. These tools are made of cast martensitic stainless steel GX30Cr13 (norm NF EN 10027-1 and Table 1). This process involves high temperature, high strain rate, thermal and mechanical stresses, fatigue phenomena, oxidation and complex moulds shape which lead to appearance of cracks on the lower mould surface. These cracks are initiated on casting defects (shrinkage).

The lifetime of the mould is determined by the dimensions of these cracks : no trace must be visible on the finished part. It as been observed that the so defined lifetime could be very different (factor 2) between two moulds of the same batch, whereas the process parameters don't seem to have vary.

This paper presents a finite element simulation of the glass pressing process to calculate the thermo-mechanical stresses in the dies. This simulation could allow to identify process parameters which small variation could involve significant stresses variation. The software FORGE2  was used because he authorizes a precise description of the cinematic of the process and possess an automatic remeshing module. To minimize time calculation, axisymetric calculation will be realized. The parameters that influence the most the stresses in the mould are : the rheological behavior of glass, the mechanical behavior of the steel, the radiative heat flux in the glass, the heat flux density at gob/mould interface and the friction between glass and mould.

Input data of the simulation

Rheological behavior of glass

The glass studied is Pyrex ® (Table 2). It is a borosilicate glass with a low expansion coefficient in order to resist up to 300°C.

Rheological behavior of glasses is very temperature dependant [START_REF] Zarzycki | Glasses and the vitreous state[END_REF]. At low temperature, glasses behave like a brittle elastic solid; around the glass transition temperature, it behave like a thermo-rheological simple visco-elastic material [START_REF] Wiertel | Optimization techniques for flat glass tempering[END_REF][START_REF] Hessenkemper | Elastic constants of glass melts above the glass transition temperature from ultrasonic and axial compression measurements[END_REF][START_REF] Van Iseghem | Residual stress simulation in a glass/metal contact[END_REF]; and at high temperature it behaves like a viscous fluid [START_REF] Lochegnies | Finite element contributions to glass manufacturing control and optimization. Part 2. Blowing, pressing and centrifuging of hollow items[END_REF][START_REF] Milutinović-Nikolić | Mathematical modelling and simulation of drawing thin glass sheet from a rectangular perform[END_REF][START_REF] Hyre | Effect of feeder design and operation on gob shape[END_REF].

Non-Newtonian behavior could be observed [START_REF] Yue | On the different descriptions of the non-Newtonian viscosity (shear-thinning effect) of glass melts with respect to heat dissipation[END_REF] : the viscosity of the glass melt decrease with the increase of strain rate. In this paper, the glass is supposed to be a Newtonian fluid (Eq.1), that means that the viscosity will be assume to vary only with the temperature as described by the Tammann-Vogel-Fulcher (TVF) equation (Eq. 2 and Fig. 2) [START_REF] Sipp | Viscosity, configurational entropy and relaxation kinetics of borosilicate melts[END_REF].

   ) ( 2 T s  (1) 0 log T T B A     (2)
Where s is the deviatoric part of the Cauchy stress tensor,  is the glass viscosity,  is the strain rate tensor, T is the temperature and A, B, T0 are constants.

Mechanical behavior of GX30Cr13 steel

It will be assumed that the steel of the mould behave as a thermo-elastic solid (Fig. 3). Young's modulus will be introduced in the software FORGE2 ® as a function of temperature (Fig. 3) and Poisson's ratio will be equal to 0.3.

Conductive and radiative heat transfer in glass

The high temperature dependency of the glass rheological behavior involves the necessity of a precise thermal computation. Glass being a semi-transparent media, heat transfer occur not only by conduction but also by radiation. Thermal conductivity (W.m -1 .K -1 ) of glass is supposed to vary with temperature (K) as described by Loulou et al. (cf. [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF] and Eq. 3).

)

273 ( 001583 . 0 7222 . 0    T k cond (3)
Radiative heat transfer is a very complex phenomena because heat flux in one point of the gob is influenced by every other point of the glass part. The easiest way to take the radiative heat transfer into account is to use a diffusion approximation : a radiative conductivity is defined as a function of the temperature, radiative and conductive heat fluxes are then calculated by Fourier's law. An expression of this radiative conductivity krad was proposed by Rosseland [START_REF] Lentes | Three-dimensional radiative heat transfer in glass cooling processes[END_REF] (Eq. 4 and Fig. 4). n is the refractive index of the semi-transparent media,  is the Stefan- Boltzmann constant, T is the temperature in K and Ross  is the Rosseland mean extinction coefficient defined by Eq. 5.  is the frequency of the radiation,

) (  is extinction coefficient and ) , ( T B 
is the blackbody intensity as described by the Planck's function. This method is very easy to compute, but the assumption is made, that the glass is optically thick. This means, that the most of the radiation is absorbed in the neighborhood of the particle that emitted it. This assumption is not valid for a glass domain with small dimensions, as in the present case.

Analytical solution of the Radiative Transfer Equation is available in the one-dimensional case, but it require to know the optical properties of the glass and the surface of the mould, and strong assumption are made on the geometry of the problem [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF][START_REF] André | A reduced analytical model for solving heat transfer in glasses[END_REF].

Numerical methods could be used as Discrete Ordinate Method, Improved Diffusion Approximation, Ray Tracing, Monte Carlo [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF][START_REF] Van Der Linden | Radiative heat transfer in hot glass melts[END_REF], but they require long calculation time.

In a first time, radiative heat transfer will be neglected because analytical and numerical methods are not easy to introduce in FORGE2 ® . Fig. 5. Comparison of profile of heat fluxes divergence in the one-dimensional case for Rosseland diffusion approximation, analytical method and without radiation Fig. 6. Heat transfer coefficient at glass/mould interface as a function of time Fig. 5 presents the profile of heat flux divergence in the glass supposed to be an infinite plate for a computed temperature profile calculated with Rosseland diffusion approximation, analytical method ( [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF][START_REF] Van Der Linden | Radiative heat transfer in hot glass melts[END_REF]) and without radiation effects. It could be observed that in this case, the Rosseland method over evaluate heat fluxes divergence compared to analytical method. In a first time, it will be assumed that radiative heat transfer could be neglected compared to conductive heat transfer.

Heat flux density at gob/mould interface

The second function of a mould being to cool the glass, the thermal computation requires to know the heat flux density exchanged between the gob and the mould. Experiments have been made [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF][START_REF] Lochegnies | Computation procedure for the temperature in hot glass. Application to the finite element simulation of hollow glass forming[END_REF][START_REF] Moreau | Experimental methodology for the determination of the heat transfer coefficient at the glass/mould interface[END_REF] to measure heat flux at the glass/mould interface as a function of time. Semi-empirical models have been developed by decoupling radiative and convective part of the flux. For example, Pchelyakov [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF][START_REF] Viskanta | Analysis of heat transfer during glass forming[END_REF][START_REF] Storck | Heat transfer modeling of the parison forming in glass manufacturing[END_REF] proposed a heat transfer coefficient equal to the conductivity divided by the thickness of a gap between gob and mould. The experimental results of Moreau et al. [START_REF] Moreau | Experimental methodology for the determination of the heat transfer coefficient at the glass/mould interface[END_REF] (Fig. 6) will be used in this paper.

Friction between glass and mould

Friction between glass and steel is largely influenced by the temperature of the glass at the interface. Studies have been led on friction between glass and steel in the condition of delivery gobs into moulds [START_REF] Falipou | New method for measuring the friction between hot viscous glass and metals[END_REF]. Li et al. [START_REF] Li | An experimental study of the lubrication behavior of A5 glass lubricant by means of the ring compression test[END_REF] investigated the influence of temperature and strain rate on the friction at the interface during compression test of titanium rings lubricated with glass. Below 950°C, strain rate does not seem to play a significant role on friction behavior. These results will be used in this paper : it is supposed that glass temperature at the interface is 750°C, and that friction behavior follows the Tresca model with a fiction factor m equal to 0.28 (Eq. 6). In Eq. 6,  is the shear stress and V   is the difference between the velocity of the glass and the velocity of the mould.

V V T m       ) (   2.

Heat capacity of Pyrex ®

Heat capacity of Pyrex ® has been introduced in FORGE2 ® as a function of the temperature has described by Richet et al. [START_REF] Richet | Configurational heat capacity and entropy of borosilicate melts[END_REF].

Simulation of the gob loading in the mould

Initial conditions

The first stage of the glass pressing process that has been chosen to simulate is loading of the gob in the mould. The initial gob shape, gob mesh and temperature distribution are presented in Fig. 7. The initial temperature distribution corresponds to a 0.5 s cooling in contact with a loading chute at 100°C on the lower side (heat transfer coefficient 2000 W.m -2 .K -1 ) and convection with air on the other ones (heat transfer coefficient 150 W.m -2 .K -1 [START_REF] Lochegnies | Inverse determination of glass viscosity and heat transfer coefficient by industrial testing[END_REF]). The gob has been placed at a distance of 50 mm from the mould and falls freely. The initial shape, mesh and temperature distribution of the mould are shown in Fig. 8. The initial temperature distribution is the result of a 30 s cooling by convection in free air (heat transfer coefficient 10 W.m -2 .K -1 ) corresponding to the time between loading the mould on the press and first pressing.

Boundary conditions

Thermal exchanges considered are convection with air during loading of the gob, convection between the mould and air, and thermal exchanges between mould and glass. The heat transfer coefficient between gob and air is equal to 150 W.m -2 .K -1 , as it can be interpolated at 1100°C from the results of the inverse method proposed by Lochegnies and Marechal [START_REF] Lochegnies | Inverse determination of glass viscosity and heat transfer coefficient by industrial testing[END_REF]. Convection coefficient between air and steel as been set to 10 W.m -2 .K -1 (the default value proposed by FORGE2 ® ). Heat transfer coefficient as been defined as a function of time (section 2.4. and Fig. 6).

Friction between glass and steel is defined in section 2.5.

Results

The total time of the simulation is 2 s. The final configuration of this simulation will be defined as initial configuration of the pressing stage (Fig. 9). Fig. 10 present the variation of the Mises equivalent stress at the surface of the mould on the axis during gob loading simulation and Fig. 11 shows the temperature at the same point. The maximum Mises equivalent stress reaches 430 MPa at 515°C (the maximum temperature reached is 520°C). Mechanical behavior of GX30Cr13 cast steel should be studied in details, but it could be supposed, that this stress could at least provoke local plastic deformation in the shrinkage zones where stress concentration may occur, and lead to cracks appearance.

Initial conditions

The second stage of the process that has been simulated is pressing. The initial configuration corresponds to the final one of the gob loading stage. Fig. 9 shows the initial temperature distribution in the gob and in the lower mould. Initial temperature of the upper mould has been set to an uniform value (250°C).

Boundary conditions

The boundary conditions for this simulation are the same as for gob loading, except the heat transfer coefficient between gob and moulds that as been set to a constant value (2000 W.m -2 .K -1 ). It has been supposed that effect of pressure on heat transfer coefficient is negligible and that his value is the final one of the gob loading stage (Fig. 6). The inner surface of the upper mould is cooled by a 20°C water flow. Fig. 12. Piloting height as a function of time.

Press piloting

The lower mould is fixed and the upper one is piloted with a height-time curve (Fig. 12). This curve has been interpolated from initial height, final height and initial velocity by an exponential function (Eq. 7). Where h is the piloted height (mm) and t is time (s).

Results

Final temperature distribution in gob and moulds is shown Fig. 13. Mises equivalent stress at the surface of the mould on the axis has been plotted as a function of time on Fig. 14 and temperature at the same point on Fig. 15. The maximum Mises equivalent stress is 850 MPa at 489°C. As it could be seen on Fig. 16, there is an other zone where stresses are important (in the edge of the plane surface of the mould). In this point, maximum Mises equivalent stress is not so important as on the axis but a high stress level exists. The evolution of the stress at this point during pressing is plotted Fig. 16. It could be concluded that as for loading stage, plastic deformation could occur during pressing around the shrinkages.

Conclusion

To validate this numerical results with experimental data, computed pressing force has been compared with the real one (the industrial process use a 5.4 tons hydraulic press). The calculated force has been plotted on Fig. 17. Although the models that have been introduced in the software are not very precise (Rosseland approximation) or corresponds to conditions that are not exactly similar to the ones of our case, the order of magnitude of the computed pressing force is realist but largely over evaluated. The difference between the calculated and the real pressing force could largely be explained by the piloting of the press in the pressing simulation which does not corresponds to measurement (the velocity is not described on a realist way). Even if the stresses in the mould are over-evaluated (as the pressing force), it could has been observed that cracks appearance could occur at two locations in the mould : at the center or at the edge of the plane surface. The first zone is solicited both thermally and mechanically, whereas the second one is almost solicited mechanically. It could has been observed that industrials moulds present cracks in one or the other of theses two zones depending on the shape of the mould : few proof and axisymetric moulds are damaged at the center whereas proof and elliptic moulds are damaged at the edge of the plane surface. This simulation could be improved by taking experimental data into account (glass rheology, mechanical behavior of GX30Cr13 steel, radiative effects, heat flux density at the gob/mould interface, piloting of the press…) and then be used to study the influence of process parameters on moulds lifetime (initial glass and mould temperature, pressing velocity…). Others mould materials could also be tried.

Fig. 1 .

 1 Fig.1. Main stages of the glass pressing process : a) delivery, b) pressing, c) pressure maintenance, d) takeout

Fig. 2 .

 2 Fig. 2. Temperature dependency of the viscosity of Pyrex ® : experimental points and VTF equation

Fig. 3 .

 3 Fig. 3. Young's modulus of GX30Cr13 steel as a function of temperature

Fig. 4 .

 4 Fig. 4. Rosseland radiative conductivity of Pyrex ® as a function of temperature

Fig. 7 .

 7 Fig. 7. Initial condition of the gob for the loading stage

Fig. 8 .Fig. 9 .

 89 Fig. 8. Initial condition of the mould for the loading stage

Fig. 11 .

 11 Fig. 11. Evolution of temperature at the surface of the mould with time during gob loading

Fig. 13 .

 13 Fig. 13. Temperature distribution in gob and moulds at the end of pressing

Fig. 17 .

 17 Fig. 17. Computed pressing force as a function of time

Table 1 .

 1

	Chemical composition (wt.%) of GX30Cr13 steel
	C	Cr	Si	Mn	Ni
	0.3	13	0.6	0.9	<0.5
	Table 2.
	Chemical composition (wt.%) of Pyrex ® glass
	SiO2 B2O3 Na2O Al2O3 K2O CaO MgO TiO2 FeO
	82.57 10.81 5.13 1.23 0.13 0.05 0.05 0.02 0.01

Acknowledgements

The present investigation has been carried out with the financial support of Aubert et Duval Alliages, Newell S.A and Saint-Gobain SEVA.