Matthieu Godichaud
email: matthieu.godichaud@laas.fr

Laas Peres

Francois ; Enit

Pérès François
email: francois.peres@enit.fr

Ayeley Tchangani
email: ayeley.tchangani@enit.fr

Optimizing end-of-life system dismantling strategy

Keywords: Uncertainty management, Bayesian networks Disassembly strategy, end-of-life modelling, Bayesian networks, optimization

The management of end-of-life systems is becoming a major concern for systems manufacturers as the negative impact of these systems on the environment is a matter of increasing public awareness and their appropriate treatment offers economic opportunities. In this context, the disassembly of these systems in order to recycle their components is a possible and sound option that can make it possible to sustain economical progress while respecting environment requirements. The work undertaken in this paper considers modelling and optimizing issues of such disassembly activities. An integrated approach is proposed to model and optimize the selection of valuable components of end-of-life systems, their recycling options and the way to obtain them. Because the framework of such problems is highly uncertain, we propose the use of Bayesian networks and their extension in terms of influence diagrams as mathematical tools for structuring and managing uncertainties. This approach allows taking into account uncertainties rising from different sources on one hand and as a support for optimization on the other hand.

 [START_REF] Presley | A strategic sustainability justification methodology for organizational decisions: a reverse logistics illustration[END_REF], the stakes related to the final step in the life cycle of systems -that is their retirement from service -have increased considerably. End-of-life systems must be disassembled in order to satisfy environmental issues. Valuable components must be selected according to technical, economic, and environmental criteria and ultimately a disassembly system has to be optimized enabling these components to be obtained.

Disassembly strategies must then respond to all decision issues raised during the stage of system retirement. In a disassembly strategy, the way to obtain and to recycle components is what we refer to as disassembly trajectory in this paper. When this trajectory is known, decision-makers can then specify the internal disassembly process, identify the needs for the repatriation logistics chain and identify and select recycling channels. A disassembly trajectory is defined with a set of products (elementary products or sub-assemblies) and a set of activities on these products. These activities could be disassembly activities, which can be destructive or not, or recycling activities which generally include:

• functional recycling that consists of introducing disassembled products into the process of new system production or into the exploitation process (maintenance for instance) of existing systems;

• material recycling that generates material from components of the end-of-life system in the production process of new systems;

• energy-oriented recycling in order to produce energy; Dangerous products and/or components that cannot be recycled may be stocked in safe places that respect environmental requirements.

Disassembly is undertaken in an uncertain context that needs to be characterized when defining a strategy. Uncertainties to be considered are mainly related to:

• the state of end-of-life systems when their deconstruction moment comes;

• the system withdrawal date;

• the withdrawal conditions and repatriation modes of the system;

• the duration of the disassembly and recycling activities and the nature of resources needed;

• the recycling channels that depend on demands of disassembled products and on their adequacy with the recycling techniques.

Our purpose is to provide the decision-makers with decision support tools that first allow the determination of disassembly trajectories based on various decision problems encountered in disassembly and second, makes it possible the management of uncertainty. Different kind of decision problems may be encountered when defining a disassembly trajectory. The first one concerns the definition of a disassembly path. It enables the determination of the intermediate sub-assemblies that have to be generated to obtain the final components. Graph theory [START_REF] Lambert | Disassembly modelling for assembly, maintenance, reuse, and recycling[END_REF], AND/OR graph and Petri nets [START_REF] Tang | Disassembly modeling, planning and application[END_REF] are the most used modeling tools to represent this kind of problem. Graph search is used to find optimal solution according to economical criteria. These models are not sufficient to represent sequence-dependant cost (i.e. cost of a sequence of operation determines the order of operation) when parallel disassembly operations need to be achieved with the same resource [START_REF] Kang | Non-linear disassembly planning for increasing the end-of-life value of products[END_REF]. This is the second kind of decision problem. It is currently modeled and solved with linear programming on the basis of a disassembly path model that specifies precedence constraints between operations (Kang et al. 2001, Lambert and[START_REF] Lambert | Disassembly modelling for assembly, maintenance, reuse, and recycling[END_REF]. In [START_REF] Moore | Petri net approach to disassembly process planning for products with complex AND/OR precedence relationships[END_REF] specific Petri nets are proposed to model this problem in which preconditions for beginning each operation are modeled. Some sub-assemblies of an end-of-life system can be recycled and do not need to be disassembled. For each sub-assembly, decisionmakers have to choose between disassembly and recycling. This leads to the third type of decision problem called disassembly depth. It necessitates the representation of the potential recycling values of sub-assemblies. Petri nets are most commonly used to represent this problem [START_REF] Zussman | A methodology for modeling and adaptive planning of disassembly processes[END_REF][START_REF] Tiwari | A Petri Net based approach to determine the disassembly strategy of a product[END_REF][START_REF] Reveliotis | Uncertainty management in optimal disassembly planning through learning-based strategies[END_REF] and AND/OR graph are also used in [START_REF] Kang | Non-linear disassembly planning for increasing the end-of-life value of products[END_REF] [START_REF] Lambert | Disassembly modelling for assembly, maintenance, reuse, and recycling[END_REF]. Indeed, these models allow integrating sub-assembly recycling values into their optimal search algorithms. The fourth type of decision problem that decision-makers may encountere in disassembly concerns the selection of a single recycling option from a set of possible options for each selected component. Petri net and associated search algorithms are also currently used to model this problem (Zussman andZhou 1999, Tiwari et al. 2002). Indeed, only a few modifications are necessary to integrate various recycling options for each component in order to adapt Petri nets used for modeling disassembly depth and path. In the disassembly of end-of-life system, destructive disassembly operation is an efficient alternative to reduce operation time and cost.

However, very few studies deal with this problem. Destructive disassembly of a component is performed when the component is too damaged [START_REF] Gungor | Disassembly sequence planning for products with defective parts in product recovery[END_REF]. In (Das and Naik 2002), the authors propose a specific process model to represent destructive operations that generate recycling components and carcasses (component or sub-assemblies that cannot be recycled). Generation of material fragments is considered in [START_REF] Kanai | Representation of product and processes for planning disassembly, shredding, and material sorting based on graphs[END_REF] as the result of destructive operations. Tree structure is used to model the material fragments that compose the system to be disassembled. In [START_REF] Duta | Determination of the Optimal Disassembly Sequence Using Decision Trees[END_REF], the authors suggest integrating destructive operations in disassembly Petri nets. Through this study, Petri nets seem to be the most flexible tool to model the different decision problems linked to the disassembly trajectory definition (Table 1).

Uncertainties are inherent in the deconstruction activity and more generally in end-oflife systems management tasks. As far as we know, few models of the literature actually take them into consideration and, when they do, very often only one type of uncertainty is dealt with. In most cases, sources of uncertainty concern the different levels of the state of degradation. In [START_REF] Krikke | On a medium term product recovery and disposal strategy for durable assembly products[END_REF], probability distribution is associated with each component and a dynamic programming procedure allows the determination of an optimal trajectory in a tree-based representation of the end-of-life system. A similar approach is proposed in [START_REF] Salomonski | On-line predictive model for disassembly process planning adaptation[END_REF], [START_REF] Reveliotis | Uncertainty management in optimal disassembly planning through learning-based strategies[END_REF] based on Petri nets and with learning approaches to determine probabilities. The same Petri net structure is used in [START_REF] Turowski | Analysis of an adaptive fuzzy system for disassembly process planning[END_REF]) but the authors model uncertainties with fuzzy logic. The failure of a disassembly operation is a consequence of product state uncertainties (degradation or constitution). A Bayesian network is proposed in [START_REF] Geiger | Probabilistic reactive disassembly planning[END_REF] to model the state of the system and the risk of failure of disassembly operations. In (Zussman and Zhou 1999, Duta et into the optimal search algorithm. Operation times and cost (disassembly and recycling) are also non deterministic since they are the result of various disassembly operations. The component states and their nature make these parameters highly uncertain [START_REF] Turowski | Analysis of an adaptive fuzzy system for disassembly process planning[END_REF][START_REF] Kang | Non-linear disassembly planning for increasing the end-of-life value of products[END_REF]). In [START_REF] Kang | Non-linear disassembly planning for increasing the end-of-life value of products[END_REF], the author also proposes to take into account the variation in recycling revenue to select an optimal disassembly trajectory. One of the principal limitations of the approaches presented so far is that they do not integrate many uncertainty sources encountered in the end-of-life systems management processes.

Table 1 sums up the related work we have analyzed. To model and solve disassembly trajectory problems, Petri nets seem to be the most flexible modeling tool. They represent every type of disassembly decision problem. Nevertheless, one of the principal limitations of the approaches based on this tool is that they do not integrate the various uncertainty sources traditionally encountered in the end-of-life system management processes. In this paper, we propose a new approach to model and determine optimal trajectory that resumes modeling flexibility of Petri nets and allows the integration of various sources of uncertainty. In section 2, a Bayesian network-based model for disassembly is presented. From this model, optimization principles are presented in section 3. An illustrative example is proposed in section 4.

Disassembly strategy modelling

Bayesian networks for disassembly modelling

It has been shown in [START_REF] Godichaud | Outils d'aide à la décision pour la sélection des filières de revalorisation des produits issus de la déconstruction des systèmes en fin de vie : application au domaine aéronautique[END_REF]) that Bayesian networks may constitute a privileged modelling tool to represent various sources of uncertainties in disassembly. In [START_REF] Godichaud | Outils d'aide à la décision pour la sélection des filières de revalorisation des produits issus de la déconstruction des systèmes en fin de vie : application au domaine aéronautique[END_REF], a modelling procedure is proposed to define a disassembly Bayesian network on the basis of a Petri net. This model seems quite interesting since it has the advantages of the disassembly Petri nets but also because it uses the capacity of Bayesian networks to deal with uncertainty management. Apart from this ability of the model to represent uncertain situations, the model makes it possible to return to the causes and advance to the consequences as well as to update probabilities based on knowledge of the context. For a precise description of Bayesian networks see for instance (Pearl 1988, Jensen and[START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF]. Let us point out however that the data required to fill in the Conditional Probability Table required to assess the situations described by the model can be obtained using statistics and / or learning as well as expert judgment.

In order to represent decision problems, decision and value nodes can be introduced in a Bayesian network. Such models are called influence diagrams. They are an extension of Bayesian networks for the representation of uncertain decision-making problems (Howard et al. 1981, Jensen and[START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF].

Three types of nodes characterize an influence diagram (see Figure 1):

• chance nodes, (generally represented by circles) represent problem variables;

• decision nodes, (generally represented by rectangles) stand for the different choices available to decision-makers,

• utility or value nodes, (generally represented by diamonds) enable the numerical evaluation of decision consequences.

Edges or arcs connecting a chance node to a decision node correspond to information available to the decision-maker when making his decision. There are various algorithms to solve decision problems represented by influence diagrams, see for instance [START_REF] Lauritzen | Representing and solving decision problems with limited information[END_REF]Nielsson 2001, Jensen andNielsen 2007).

Problem representation

The disassembly trajectory problem is represented by Bayesian networks (BN). Indeed, they enable all the elements of this decision problem to be represented. Generally speaking, system disassembly modelling with Bayesian networks is described by the following items: • "product" nodes representing end-of-life system components that have one or more recycling option

• "activity" nodes representing disassembly operations or recycling action on each product,

• arcs characterizing precedence and exclusion relationships between activities,

• node parameters that make it possible to characterize disassembly process progress. Based on the generic representation of an activity in a process, an example of disassembly trajectory is presented on Figure 2. The deconstruction of a system modeled by node S C is considered. This node corresponds to the input flow of the activity modeled by node A 1 . It is a disassembly operation which generates the products modeled by nodes P 2 and P 3 . Product P 2 can be recycled by carrying out the activity A 2 . Two activities are then possible on product P 3 : material recycling (node A 3 ) or functional recycling (node A 4 ). Only one activity must however be selected. A decision node C 1 represents the selection between both activities.

Purposes of the disassembly model are:

• representing the disassembly processes,

• taking into account various origin uncertainties,

• evaluating disassembly trajectories.

Modelling of the disassembly process

The modality set of an "activity" node (disassembly operation and recycling action) corresponds to the various realization modes of the modeled activity. In the most straightforward configuration corresponding to the description of the chaining of the process activity, two modalities are necessary:

• " ": activity is carried out,

• " ": activity is not carried out.

Main input and output elements in a disassembly trajectory correspond to the components of an end-of-life system (sub-assemblies or elementary components). To describe the realization of a disassembly trajectory, a minimum of two modalities is necessary for the "product" node corresponding to an input or output element: • " ": the element is activated i.e. the beginning conditions of an activity are realized,

• " ": the element is not activated.

The basic mechanism to be modeled is characterized by the repetition of the following phases: product activation and activity realization. Product activation is modeled by the arcs going from "activity" nodes to "product" nodes. The mechanism to be modelled works as follows:

• when the activity is not realized (modality " "), the output products are not generated or activated (modality " "),

• when the activity is carried out (modality " "), the output products are generated or activated (modality " ").

In cases where there is no uncertainty, the element activation is modelled in conditional probability tables (CPT) of the "product" nodes as illustrated in Table 2 which must be read in the following way:

• , i.e. the "product" node takes the modality " " with probability equal to 1 when the activity node takes modality " " (deterministic case),

• , i.e. the "activity" node takes the modality " " with probability equal to 1 when the activity node takes the modality " ".

Following the same principles, activity realization is modelled by "activity" node CPT according to arcs going from "product" nodes to "activity" nodes. An activity can be realized only if the product to be transformed is activated. The specification of the activity realization modelling is presented in Table 3.

Handling disassembly uncertainties

System state uncertainties

System state uncertainties are taken into account in "product" nodes. Degradation levels of each component of an end-of-life system have to be characterized. To consider various possible states of the products obtained through disassembly operations, the set of "product" node modalities are modified. The modality " " (not activated) is used to model the process realization progress. Modality " " is modified for products that can have different degradation levels. A "product" node will be thus characterized by the modality set if degration levels have to be taken into account.

Disassembly operation uncertainties

Another uncertainty source is related to the nature of disassembly operation. Systems are not necessarily designed to be disassembled and operations are not standardized. Thus there are intrinsic uncertainties when determining disassembly effort. To handle uncertainties relating to disassembly operation realization modes, the modality set of "activity" nodes is modified.

Modality " " is used to model the process realization logic. Modality " " can be replaced • operation failure when the operation duration becomes too important or when it is not realizable ;

• : not selected operation.
There is no generic situation for parameter specifications on "activity" nodes. Nodes representing disassembly operations can indeed be conditioned by various other nodes according to the operation environment. Three current situations can however be highlighted:

• direct evaluation: uncertainty is relating to not modeled factors i.e. the node is not conditioned by another: the user will evaluate the probabilities ,

• evaluation according to product conditions: state of disassembled product can influence operation realization ; probabilities associated with a node are then conditioned by a node : the user has to evaluate ,

• evaluation according to the resource used: when various resources can be used for an operation, they can influence its realization: if represents a resource, the user has to evaluate , where is the number of possible using mode of the resource.

Recycling action uncertainties

Recycling action uncertainties are taken into account on the associated "activity" nodes and, if necessary, on nodes representing the demand for disassembled products. Modelling principles of increasing value actions are the same as those used for the disassembly operations. An Modelled uncertainties can, for instance, being related to demand dates or cancellation. In this case, a node modelling the demand for a product will have modality set where corresponds to possible demands: corresponds to demand cancellation and may represent the probabilities of demand at a given period.

Economic evaluation

Disassembly strategies can be evaluated from an economic point of view through elements relating to:

• disassembly operation costs,

• recycling action costs (realization and storage),

• recycling value incomes.

These economic elements are modelled by utility nodes attached to each disassembly process activity as illustrated in Figure 3 . We can retrieve the utility form of Figure 7(a)

by applying the equation [1] (used in optimization described in the following section):

The purpose of the model proposed in this section is to represent the problem of disassembly trajectory determination i.e. a set of activities with their connexions. It also enables these various activities to be analyzed and specified in order to serve as a support for optimization mechanisms which are presented in the next section. The optimisation mechanism consists then in tracking each product model in a recursive way going from the elementary component models to the global end-of-life system.

Product disassembly generic model

The model presented above represents the whole disassembly trajectory that the decision maker has identified. The model is a network and the objective is to find an optimal trajectory within this network that for each product, given its state, what the best activity would be. In this network, represents the set of successors of "product" node (i.e. )

and represents the set of product node successors of a disassembly activity . A disassembly policy model is drawn from the global model to evaluate each product separately.

It enables the required defining recursive equation to be obtained to determine the optimal disassembly trajectory.

Disassembly policies are modelled by decision nodes associated with each product.

These nodes are integrated in the model as presented in Figure 4 (node ) which gives a generic model representation (integration of all the elements required to determine a policy).

The considered product is modelled by node and modalities of node characterize all the possible options likely to be selected on the product. Utilities (i.e. is a component of ) represent the evaluation of product components generated by each disassembly operation. A policy model being associated with each product, these utilities correspond to the optimisation result of the product component policies.

Reduction rules for optimization

Given the modelling generic structure proposed in Figure 4 and the model specificity (activation principles of variables), the method using diagram progressive reduction is appropriate. Indeed it shows how recursive techniques can be applied to evaluate trajectories by considering each valuable component of the system. Consequently, resolving and analyzing the disassembly of a given system is made easier. We introduce the resolution technique in this section in order to apply it to the generic model used for the assessment of a product valuation policy.

The resolution method is made up of four generic operations (Jensen and Nielsen,07) that are applied to the influence diagram. These operations enable the progressive reduction of the model and step by step resolution until only one utility node remains. This node stands for the expected utility of the end-of-life system, given that at each decision node the best solution has been selected. The four reduction operations are: merging utility nodes, removing chance nodes, removing decision nodes, changing arc direction. The use of these four rules enables the resolution of any influence diagram:

• (R1) Merging several utility nodes: Let us consider a set of utility nodes … of an influence diagram. These nodes can be merged into a single utility node U * whose conditional expectation is given by :

• (R2) Removing chance node: Let us consider a node father of a single utility node . Given the node is represented by the probability distribution and node is described by its conditional expectation , node can be removed by modifying the conditional expectation of node as followed:

• (R3) Removing decision nodes: Let us consider a decision node father of a single utility node whose predecessors are also predecessors of node (predecessors of node are presumed to be known at the time the decision is made). Given that the policy relating to node has to be evaluated and that node U is described by its strategy is then defined as follows

• (R4) Changing arc direction: when none of these rules can be applied, some arcs can be reversed in order to reach to a situation where the previous rules can be used. To this end, rules belonging to the probability theory must be applied. Some examples are presented in (Tatman andShachter 1990, Jensen andNielsen 2007).

Application to the disassembly model

Given the policy model presented in Figure 4, evaluation of a product in state and its optimal policy are defined by equations [6] and [7], with :

If the costs and/or incomes of an activity are dependent on product state, utility must be replaced by . If the costs of a recycling action are dependent of a demand for the product then equation [1] could be used to retrieve the shape of equations [6] and [7]. The generic resolution method is made of different steps corresponding to the application of a reduction rule on a variable or a set of variables of same nature. These different steps are the following:

(S1) Removing nodes representing the recycling action;

(S2) Removing utility nodes linked to nodes representing the disassembly operations;

(S3) Removing nodes representing the dismantling operation;

(S4) Merging utility nodes relating to the given policy;

(S5) Removing decision node representing the policy;

(S6) Removing node representing the product to be recycled.

Input data correspond to:

• probabilities for all activities (disassembly and recycling) that can be performed on ,

• probabilities for all components of generated by disassembly operation • utilities values for disassembly operation that correspond to disassembly operation costs,

• utilities values for recycling activities that correspond to difference between incomes and recycling action cost.

The first step (S1) corresponds to the determination of the utility of a recycling action given that product is in state . These nodes can be removed since they have as single successor, utility nodes (rule R 2 ). Parameters are modified as follows:

If is an elementary component, it is associated with recycling action only. Utilities and the corresponding policy can be determined directly after this step.

The second step (S2) relates to the determination of the utility of a disassembly operation given that product is in state . It depends on disassembly cost and utilities for all components generated by . Utility nodes and connected to the same node representing a disassembly operation can then be merged (rule R 1 ):

At this step, nodes representing disassembly operations can be removed, given that they have as successors only utility nodes (rule R 2 ). Utilities of each disassembly operation on is then given by (step (S3)): The set of utility nodes has now the same parents and . They can be merged in the fourth step into a single utility node (rule R 1 ). The fifth step (S5) consists of determining the activity for P. Utility and policy for P is then:

It corresponds to [6] and [7] if we express with equation [8], [9]
and [10]. Eventually, (S6) corresponds to the removal of node .

Illustrative example

System layout

The problem we want to tackle concerns the valuation trajectory of an airplane turbine represented in Figure 6. The system is made up of seven elementary products (P1 to P7) and sub-assemblies are: SA1 = {P2, P3, P4, P5, P6} ; SA2 = {P2, P4, P5} ; SA3 = {P3, P5, P6}. The connection graph on Figure 6(b) is used to identify the disassembly operation and the precedence relationships. The method of analysing a connection graph is presented in detail in [START_REF] Lambert | Disassembly modelling for assembly, maintenance, reuse, and recycling[END_REF] for instance and its application to the example is developed in [START_REF] Godichaud | Outils d'aide à la décision pour la sélection des filières de revalorisation des produits issus de la déconstruction des systèmes en fin de vie : application au domaine aéronautique[END_REF]).

For the system considered here, all the elements are recapitulated in Table 5 which gives a rapid overview of the model nodes. Disassembly operations concern only subassembly and two types of recycling action are considered: material recycling and functional recycling. The states of the nodes of the model are defined in Table 6 in a generic way for each type of node. The states of the decision nodes characterize the possible activities for each product. Consequently, a decision node representing a policy has a number of states standing for the different activities associated with the considered product.

Utility nodes stand for the performances which depend on a costing model, itself based on the parent nodes. To each configuration of the parent nodes a value is defined.

The model input parameters correspond, on the one hand, to the uncertainties relating to activity realization as well as to the state of the products and, on the other hand, to activity costs and incomes. All the model input data of the valuation trajectory is presented in Table 7 for the "activity" and "utility" nodes. Table 8 shows the "product" nodes.

Optimization result

An algorithm based on the recursive relations [6] and [7] has been implemented with the graph of Figure 7 and its parameters as inputs. It proceeds by beginning the evaluation with elementary components P1 to P7 and the sub-assembly SA2 which cannot be disassembled.

Sub-assembly SA3 can then be evaluated since all its components were evaluated. It is then the sub-assembly SA1 which can be evaluated and finally the complete system.

The graphical representation of the optimal trajectory is given on Figure 7 products P 1 , SA 1 , P 7 , P 2 , P 4 and SA 3 . For each product, the optimal strategy is indicated in the decision nodes (rectangles). For product SA 1 , for instance, it shows that activity DO 3 has to be selected. For product P 1 and P 2 , it depends on their state. If the products are in a deteriorated state "d", the trajectory is indicated with black non-bold lines (the non-selected trajectory being represented by a dotted line).

The result of optimisation is presented on The optimal trajectory is presented on Figure 8 with the probabilities (percentages) at each node. For each "product" node, the values correspond to the probabilities of generating this product in various states (normal or degraded) and the probabilities of not generating this product (especially when the trajectory fails upstream). For each "activity" node, the values correspond to the probabilities of realizing this activity in various modes (normal or slow), stopping it or not realizing it. For instance, if 100 arrivals of systems S are expected on a given period with a probability of 0.5 of having a normal state or degraded state, the product P 7 is expected to be 92.4 times in a normal state, 5.12 times in a degraded state and it is not generated 2.5 times. Furthermore, the recycling activity FR 9 of P 7 is expected to be realized 67.7 times in a normal mode, 24.6 times in a degraded mode and stopped 5.13 times (not realized 2.5 times).

Conclusion

The end-of-life system disassembly strategies are based on the determination of trajectories completely defined by the products subject to increasing value actions, their mode of obtaining as well as their valuation circuits. Trajectory modelling and optimizing methods within an uncertain context have been proposed in this article. The model structure provides the guidelines for analyzing the trajectories. It also enables the consideration of uncertainties of various natures in the problem representation and optimization as they generally appear in the dismantling context of end-of-line systems.

The working perspectives concern first the consideration of uncertainties likely to evolve with time (duration of an activity, date of arrival or request …). This type of modeling will enable the planning of several arrivals of end-of-life systems to be disassembled. Within this framework, Bayesian dynamic networks provide an interesting modeling solution.

We have used an economic criterion to evaluate the quality of a disassembly trajectory. Other criteria must be taken account in order to have better control over the disassembly process. These are environmental criteria based on the ecological impact of the disassembly process. Multi-criteria optimizing approaches could be used for this purpose.

Another perspective deals with the model deployment on a real industrial case. At the moment the model is being applied in the field of airplane recycling.

  of sustainable development concerns or considerations

  ), a probability of failure is associated with each disassembly operation and integrated

  Decision variables are attached to each product. They indicate the direction of the disassembly trajectory towards one option (disassembling or recycling). Constraints are specified by the arcs. Economical parameters are associated with "activity" nodes by means of utility nodes. They represent costs and incomes potentially generated by the realization of an activity. They enable the economic profit of the various trajectories to be evaluated. The set of nodes of the Bayesian networks disassembly model is noted . The following subsets of nodes characterize the model structure: • is the set of "product" nodes with : o an element of , o the node representing the whole system o the set of modalities of node , • is the of "activity" nodes with : o subset of "activity" nodes representing disassembly operations and representing disassembly operations on , o subset of "activity" nodes representing recycling action and representing recycling action on , o an element of , o the set of modalities of node , of utility nodes. They are associated with each activity and is the utility node associated with the activity modelled by node .

  is associated with each identified action. Its modality set represents the activity realization modes. These various modes characterize the duration as well as the required resources needed to recycle the product.Recycling action generates an income if the related product is subject to a demand. Demands are determined by forecasts which can be uncertain. Nodes characterizing the demand are then integrated into the disassembly model. Demand uncertainties imply considering storage of product inducing handling costs as presented in the next section.

  These equations can be justified by a step-by-step analysis of the policy model. Let us apply the different reduction rules previously introduced (equations [2] to[5]) in order to optimize the dismantling trajectories represented by influence diagrams. To this end we will proceed to the optimization of the generic model of Figure4. This model can be used at each decomposition level of the end-of-life system. The evaluation method consists of processing each elementary model in a recursive way by going from the basic components to the product to be recycled as presented on Figure5. The generic resolution method is made up of different stages corresponding to the application of a reduction rule on a variable or a set of variables of the same nature.

  the model Nodes are interconnected according to the system structure and the relation highlighted in the connection diagram in Figure6. The general form of the model is presented in Figure7. Each identified product has the same generic structure introduced in Figure4. It is customized according to the disassembly operations and recycling action likely to be performed on the product. For the SA 1 product for instance, one comes across the nodes corresponding to the recycling action MR 1 and FR 1 as well as the nodes relating to the disassembly operations DO 2 and DO 3 .

  in bold. The optimal trajectory is drawn in bold lines. It shows that the best solution consists of obtaining

FigureFigure 6 .Figure 5 .Figure 7 .

 657 Figure 4. Disassembly policy model

  . Economic parameters are then specified in the utility table attached to each utility node.

	activity modelled by	according to its different realization modes. The general shape of a
	utility table	is given in Table 4.
	The model enables other elements that have an impact on costs to be taken into
	account, such as product state conditions or a particular context (demand for the product for
	instance). Different situations are presented in Figure 3(b) and (c). Utility table elements
	correspond realization modes and to the state of the product or the demand respectively. If then to and F , namely the evaluation of activity according to activity , with o r ) is
	P the realization cost of activity	in mode	and	is the income generated by recycling
	e through , utility values associated with e r incomes depend on product will be	state,	the	utility	. If recycling values
	become			R	. In Figure 3(c), demand has an impact on
	e e i v recycling cost because there are storage costs. If	) is the storage cost then
						w
						O n l
						y
	A utility function is associated with each utility node. This function gives a value for
	each configuration of parent nodes of the considered utility node. If a node	is considered
	with the set	of its parent nodes, utility of	is noted	. A utility	of a
	utility node associated with an "activity" node , characterizes costs and incomes of the

  Table9for each product. All the products (component and sub-assemblies) are listed in the first column. The second column states whether or not each product is in the optimal disassembly trajectory. The third and fourth columns give the policies for each product given their states (normal state PL(P=n) or degraded state PL(p=d)). For instance, the policy for S consists in realizing disassembly operation DO 1 whatever its state. Utility values for each product given its state are presented in the last column. The disassembly utility of S is 243.3 (monetary unit here because costs have been considered as utility values) if it is in normal state and 231.23 otherwise. If 100 arrivals of systems S are expected at a given period with a probability of 0.5 of having a normal state or degraded state then the expected incomes for this period are 23758.

Table 2 .

 2 CPT of a "product" node with an "activity" node as an input.

									S
								n	d	na
								50	50	0
									DO1
								n	s	g	nr
								92.5	5	2.5	0
			P1 d 91.9 5.62 2.5 n na		SA1 d 92.6 4.88 2.5 n na			P7 d 92.4 5.12 2.5 n na
		MR3			FR3			DO3			FR9
	n	s	g	nr	n	s	g	nr	n	s	g	nr	n	s	g	nr
	3.94	1.41	0.28	94.4	68.9	13.8	9.19	8.12	87.3 5.12 5.12 2.5	67.7 24.6 5.13 2.5
					P2					P4				SA3
				n	d	na			n	d	na			n	d	na
				82.4 10 7.62		73.4 19 7.62			87.5 4.88 7.62
		MR4				FR4			FR6			FR2
	n	s	g	nr	n	s	g	nr	n	s	g	nr	n	s	g	nr
	8	2	0	90	49.4 16.5 16.5 17.6	44.3 44.3 3.8 7.62	64.2 14.1 14.1 7.62
	Figure 8. Optimal trajectory probabilities

Table 3 .

 3 CPT of an "activity" node with a "product" node as an input.

	...	...
	...	...	0

Table 4 .

 4 Utility table: general shape.

	Product	Disassembly operation	Material recycling	Functional recycling
	S	DO1	-	-
	SA1	DO2 ; DO3	-	FR1
	SA2	-	MR2	-
	SA3	DO4	-	FR2
	P1	-	MR3	FR3
	P2	-	MR4	FR4
	P3	-	MR5	FR5
	P4	-	MR6	FR6
	P5	-	MR7	FR7
	P6	-	MR8	FR8
	P7	-	MR9	FR9

Table 5 .

 5 Nodes used in the model.

	Nodes	Meaning	Notations
			n : normal
	Product	Deterioration level	d : deteriorated
			na : not activated
			n : nominal
	Activity	Activity realization mode	s : slow g: giving up
			nr : not realized

Table 6 .

 6 Possible node values.

		Pr(A=n/P=n) Pr(A=s/P=n) Pr(A=g/P=n) Pr(A=n/P=d) Pr(A=s/P=d) Pr(A=g/P=d) U(A=n) U(A=s) U(A=g)
	DO1	0.95	0.05	0	0.9	0.05	0.05	-5	-7	-10
	DO2	0.9	0.1	0	0.8	0.1	0.1	-5	-10	-15
	DO3	0.9	0.05	0.05	0.8	0.1	0.1	-8	-10	-12
	DO4	0.95	0.05	0	0.9	0.05	0.05	-10	-15	-20
	MR1	0.3	0.3	0.4	0.3	0.3	0.4	30	10	-5
	MR2	0.9	0.1	0	0.8	0.2	0	60	30	-3
	MR3	0.8	0.15	0.05	0.7	0.25	0.05	35	15	-5
	MR4	0.8	0.2	0	0.8	0.2	0	20	10	-2
	MR5	0.8	0.1	0.1	0.7	0.2	0.1	30	15	-2
	MR6	0.9	0.1	0	0.8	0.1	0.1	20	10	-1
	MR7	0.9	0.1	0	0.8	0.15	0.05	20	10	0
	MR8	0.9	0.1	0	0.8	0.1	0.1	20	10	-2
	MR9	0.8	0.15	0.05	0.75	0.15	0.1	30	10	-5
	FR1	0.6	0.2	0.2	0.2	0.3	0.5	50	25	-10
	FR2 FR3 FR4 FR5 FR6 FR7 FR8	0.7 0.75 0.6 0.9 0.5 0.8 0.6	0.15 0.15 F 0.2 0.1 0.5 o 0.1 r 0.4	0.15 0.1 0.2 0 0 0.1 0	0.6 0.2 0.4 0.5 0.4 0.5 0.5	0.2 0.3 0.3 0.25 0.4 0.25 0.25	0.2 0.5 0.3 0.25 0.2 0.25 0.25	200 40 30 40 60 45 50	35 20 10 30 10 35 25	-4 -3 -5 -6 -10 -4 -5
	FR9	0.7	0.25	0.05	0.6	0.3	0.3	50	20	-5
			P
			e
				e r
					R
					e
					v i e
						w
						O n l
								y

Table 7 .

 7 "Activity" and "utility" node parameters.

			Pr(P=n/A=n) Pr(P=d/A=n) Pr(P=n/A=s) Pr(P=d/A=s)
		P1	0.95	0.05	0.8	0.2
	DO1	SA1	0.95	0.05	0.95	0.05
		P7	0.95	0.05	0.9	0.1
		SA2	0.9	0.1	0.85	0.15
	DO2	P6	0.9	0.1	0.9	0.1
		P3	0.9	0.1	0.85	0.15
		P2	0.9	0.1	0.75	0.25
	DO3	P4	0.8	0.2	0.7	0.3
		SA3	0.95	0.05	0.9	0.1
		P6	0.95	0.05	0.9	0.1
	DO4	P3	0.95	0.05	0.8	0.2
		P5	0.95	0.05	0.9	0.1
	Table 8. "Product" node parameters.
	P		PL(P=n) PL(P=d) U(P=n) U(P=d)
	S		DO 1	DO 1	243.93	231.23
	SA 1		DO 3	DO 3	177.63	167.44
	SA 2		MR 2	MR 2	57	54
	SA 3		FR 2	FR 2	144.65	126.2
	P 1		FR 3	MR 3	37.2	28
	P 2		FR 4	MR 4	19	18
	P 3		FR 5	FR 5	39	26
	P 4		FR 6	FR 6	35	26
	P 5		FR 7	FR 7	39.1	30.25
	P 6		FR 8	FR 8	40	20
	P 7		FR 9	FR 9	39.75	34.5

Table 9 .

 9 Product evaluation results

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research


Rana Elias

Nelly Point
email: nelly.point@enpc.fr

Joanna Bodgi

Pierre Argoul

Vibrations, Shocks and Noise How To Retrieve The Normal Modes Using The POD

Keywords: dynamic system, POD, normal modes, accuracy of an approximation

Proper orthogonal decomposition (POD) is an increasingly popular way to analyze data, and to obtain either a lowdimensional approximated description of a high-dimensional process, or useful information for damage assessment. In the case of a dynamic system with n degrees of freedom, the purpose of POD is to retrieve the modal properties from the measured response. Until now, POD has been used for systems with a diagonal mass matrix. The aim of this presentation is double: first, to demonstrate that POD can also be used for a non-diagonal symmetric mass matrix; and second, to present sufficient conditions on the response sampling, in order to retrieve the modal characteristics with a prescribed accuracy. The conditions to obtain this approximation with a given accuracy are first explicitly given for the case without damping. Then the case of proportional damping is addressed and similar conditions are shown. The critical case of two modal frequencies close together is also studied, given that it requires particular conditions. The obtained conditions show that the expected accuracy is explicitly limited by the damping ratio. Some numerical tests illustrate the accuracy evolution of the approximated normal modes obtained by this method, with respect to the variation of the observation time and the damping ratio. This careful analysis can be useful for finding the cause of poor approximation properties in more complex cases, such as analysis of variation in nonlinear normal modes.

Introduction

Studying the behavior of structures in various common situations such as earthquakes, high wind, etc., is very important for civil engineering. These studies are based on the mode shapes and modal frequencies of the structures, explaining why the development of new efficient mathematical techniques to retrieve this information is so important in the scientific field.

In this paper, we study a mathematically based technique that allows one to find the mode shapes and modal frequencies of a structure depending on the type of initial data. This method uses Proper Orthogonal Decomposition (POD) in order to determine the modal characteristics of a second-order free linear discrete dynamic system with n degrees of freedom (DoF). While several previous works have proven the accuracy of this method when the mass matrix is proportional to the identity matrix, in this paper we will prove its efficiency for all symmetric, definite and positive mass matrices, and even in the case of light proportional damping.

Two problem types will be studied in this paper: the direct problem and the indirect problem.

In the direct problem, the mass matrix M, the damping matrix C and the stiffness matrix K are known, and are used to determine the displacement vector x, the eigenvectors and eigenvalues of the system and their properties.

In the indirect problem, the displacement matrix X is used to retrieve the modal characteristics of a structure. Engineers very often have to deal with this type of problem. Indeed, when working on real structures, the matrices M, C and K remain unknown. Determining the modal characteristics of the system thus cannot be completed with the methods used for a direct problem. Therefore, it is valuable to develop methods that use only the displacement matrix X, which can be easily obtained through experimental data.

To find the best approximation of a system's characteristics, one relies on the POD: the eigenvectors of the correlation matrix R approach the eigenmodes of the system when M is a scalar matrix. When using this method, scientists found severe limitations of the POD. The first problem is that the eigenmodes of the system are orthogonal with respect to the mass matrix M, whereas the eigenvectors of the system are orthogonal to each other. In fact, the approximations of the eigenmodes are the eigenvectors of the matrix R * M. Studies in [START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF] and [START_REF] Golinval | Physical Interpretation of The Proper Orthogonal Modes Using the Singular Value Decomposition[END_REF] showed that when the mass matrix of a structure is proportional to the identity matrix, the eigenvectors of R converge to the modes of vibration. The accuracy of the method has been verified by a large number of numerical tests. It has been proven numerically that, using the matrix R * M, the eigenmodes can be retrieved even when the mass matrix M is not proportional to the identity, but to a known diagonal matrix [START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF]. A second problem appears when two frequencies are close to each other [1(p.1-7)], since POD cannot distinguish between them.

The purpose of this work is to extend the use of the POD method in order to retrieve the modal characteristics of a structure for indirect problems, in the case when the mass matrix is symmetric and positive definite, but not necessarily diagonal. Moreover, quantitative conditions on the total observation time are given in order to approximate the eigenmodes with a given accuracy. It is also shown how these conditions are modified when the system has frequencies that are very close together. These conditions are also given for a low proportionally damped system.

Researchers have recently developed this method for different domains. POD was successfully adapted for discrete linear symmetric undamped systems [START_REF] Golinval | Physical Interpretation of The Proper Orthogonal Modes Using the Singular Value Decomposition[END_REF] and for light damping systems [START_REF] Feeny | Interpreting Proper Orthogonal Modes of Randomly Excited Linear Vibration Systems[END_REF]. This approach also became useful for the case of nonlinear systems [START_REF] Bergman | The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems: An Overview[END_REF][START_REF] Kerschen | On the Model Validation in Non-Linear Structural Dynamics[END_REF], discretely sampled continuous systems (uniform sampling), [START_REF] Feeny | On the Proper Orthogonal Modes of Continuous Vibration Systems[END_REF] and even (with a modified version of POD) for a non-uniform discretization with homogenous structures.

In the literature, one can find three closely related methods that can be considered to be POD methods: the Karhunène-Loève decomposition (KLD) (for continuous time), the principal component analysis (PCA) (for random discrete variables), and the singular value decomposition (SVD) (which is originally a mere algebraic decomposition) [START_REF] Lee | Proper Orthogonal Decomposition and its Applications -Part I: Theory[END_REF].

This presentation is divided into five main sections. In section two, the direct problem is studied. It is shown how one can retrieve the modal characteristics of a structure in a rather simple way. In section three, we are interested in solving indirect problems. A POD-based method is proven to be accurate for all symmetric positive definite mass matrices. Conditions on the sampling frequency and the total observation time are developed in order to approximate the eigenmodes and the eigenvectors of the system, with a given accuracy ε. Both damped and undamped systems are studied. Numerical results are presented in section four, and a conclusion is developed in the last section.

Direct problem

This dynamic system considered here is governed by the following second-order differential equation:

(t)

+ (t)+ (t)= M x C x Kx 0 ɺɺ ɺ (1)
with 0 ( 0 ) = x x and 0 ( 0 ) = x x ɺ ɺ and with M the mass matrix, C the damping matrix and K the stiffness matrix. The matrices M, C and K are symmetric and positive definite. The purpose of this section is to determine the eigenmodes φ φ φ φ of the system. Therefore, the system given by ( 1) is transformed into a system of n uncoupled ordinary equations [START_REF] Lee | Proper Orthogonal Decomposition and its Applications -Part I: Theory[END_REF]. Applying the coordinate transformation, as shown in [START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF][START_REF] Feeny | On the Proper Orthogonal Modes of Continuous Vibration Systems[END_REF] and in [5(p.67)] :

1 2 ( ) ( ) t t - = x M q (2)
to the equation (1) yields to:

1 1 1 2 2 2 ( ) ( ) ( ) t t t - - + + = .. . M q CM q KM q 0 (3)
Multiplying by 1 2 -M recasts the system as:

1 1 1 1 2 2 2 2 ( ) ( ) ( ) t t t - - - - + + = .. . q M CM q M KM q 0 (4)
Since the matrices M, K and C are symmetric and positive definite, the matrices It should be noted that a symmetric matrix is always diagonalizable over an orthogonal basis of eigenvectors. For proportional (Rayleigh) damping, the matrices

M M C

have the same eigenvectors i p , so these two matrices can be diagonalized by the same change of basis. Only the proportionally damped case is considered in the following. Let P be the matrix the columns of which are the orthogonal vectors i p . By setting, ( ) ( )

t t = q P r (5)
and multiplying by T P , the system yields:

( ) ( ) ( ) ( ) ( ) i i t t t + + = r d i a g r d i a g r 0 ɺɺ ɺ µ λ (6)
where i µ and i λ are respectively the eigenvalues of

1 1 2 2 - - M M C and 1 1 2 2 - - M KM . Since 1 1 2 2 - - M KM
i i i i i λ ω µ ω ζ = =
Thus, the system (6) becomes:

2 ( ) (2 ) ( ) ( ) ( ) 0 i i i t t t ω ζ ω + + = r d ia g r d ia g r ɺɺ ɺ (7)
Since this above system is completely uncoupled, it is possible to solve an ordinary differential equation for each {1, ..., } i n ∈

. The characteristic equation for ( ) i r t is given by: 2 2

( 2 ) 0

i i i s s ω ζ ω + + = (8)
If the damping is low, i.e., 1 i ζ < , the equation ( 8) admits two complex roots defined as follows:

2 1 i i i i s i ω ζ ω ζ = - ± - (9)
By setting

2 1 i i i ω ω ζ = - ɶ
, the system response due to initial conditions may be expressed as:

[ ] 1 ( ) (0) cos( ) (0) (0) sin( ) i i t i i i i i i i i i r t e r t r r t ω ζ ω ω ζ ω ω -   = + +     ɶ ɶ ɺ ɶ (10)
Let i a be defined by: ( )

2 2 2 2 1 (0) (0) (0) i i i i i i i a r r r ω ζ ω   = + +     ɺ ɶ (11) In the undamped case, i a is twice the energy 2 i σ .

By introducing the angle 0

i θ π ≤ < , i.e., [ ]
1 (0) (0) tan (0) i i i i i i i r r r ω ζ ω θ - + = ɺ ɶ (12)
one can express ( ) i r t as follows:

( ) cos( )

i i t i i i i r t e a t ω ζ ω θ - = - ɶ (13) Since 1 1 2 2 ( ) ( ) ( ) t t t - - = = x M q M Pr , introducing φ φ φ φ , such that: 1 2 M P - -- - φ = φ = φ = φ = (14)
the response ( t )

x becomes:

n i i i 1 (t) r(t) r (t) x = = = = = φ = ϕ = φ = ϕ = φ = ϕ = φ = ϕ ∑ ∑ ∑ ∑ (15)
where the i φ are the columns of φ φ φ φ . These are the eigenmodes of the system. Therefore, as

M M C

(case of proportional damping), then:

1 1 2 1 2 2 2 i i i i i i - - - ⇔ = ω = ω M KM p p M Kφ φ and, 1 1 1 2 2 i i i i i i i i 2 2 - - - = ⇔ ω ζ = ω ζ M CM p p M Cφ φ
This emphasizes that ( i φ ) form a basis of eigenvectors for

T T 2 T i i i , = ( ) an d = (2 ) M I K d ia g C d ia g φ φ = φ φ ω φ φ ω ζ φ φ = φ φ ω φ φ ω ζ φ φ = φ φ ω φ φ ω ζ φ φ = φ φ ω φ φ ω ζ
In addition, from ( 14), one has

1 2 T T P M - -- - φ = φ = φ = φ =
, and so:

1 1 2 T T P M M - -- - φ = = φ φ = = φ φ = = φ φ = = φ (16)

Indirect problem

In practice, the mass, stiffness and damping matrices are unknown, the only given information about the structure is the displacements (eventually the accelerations), and the data are usually discretized in time. The main purpose of the POD is to find the eigenmodes of the initial differential system by finding the eigenmodes of a matrix R called the correlation matrix.

For m samples in time of the n-dimensional vector x, the (n x m) data matrix X is given by:

1 1 1 m n 1 n m x (t ) x (t ) x (t ) x (t ) X ⋮ ⋱ ⋮ … … … …                 = = = =                 … … … …         (17)
and the correlation matrix is defined by:

1 m T R XX = = = =
As seen in equation (15):

n i i i 1 (t) r(t) r (t) x = = = = = φ = ϕ = φ = ϕ = φ = ϕ = φ = ϕ ∑ ∑ ∑ ∑ hence, n k i k i i 1 (t ) r (t ) x = = = = = ϕ = ϕ = ϕ = ϕ ∑ ∑ ∑ ∑
, and therefore, the data matrix X can be written as:

1 1 1 m n 1 n m r (t ) r (t ) r (t ) r (t ) X ⋮ ⋱ ⋮ … … … …                 = φ = φ = φ = φ                 … … … …         (18)
Consequently,

1 m T R XX = = = = is expressed by: ij 1 r m T R         = φ φ = φ φ = φ φ = φ φ         (19) with 1 1 1 m 1 1 n 1 ij n 1 n m 1 m n m r (t ) r (t ) r (t ) r (t ) r r (t ) r (t ) r (t ) r (t ) ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ … … … … … … … …                                         = = = =                                         … … … … … … … …                 (20)
where, for 1 i, j n

≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤ , the general term i j
r of the i-th row and j-th column has the following expression :

1 1 ( ) ( ) m i j i k j k k r r t r t m = = = = = = = = ∑ ∑ ∑ ∑ (21)
When the mass matrix M is proportional to the identity, numerous studies [START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF][START_REF] Golinval | Physical Interpretation of The Proper Orthogonal Modes Using the Singular Value Decomposition[END_REF] have numerically shown that the eigenvectors of the correlation matrix

1 m T R XX = = = =
are a good approximation of the system's eigenvectors. It will be proved in the following section that this result is also verified even when M is no longer proportional to the identity and not even necessarily diagonal.

The main step is to show that under certain conditions, the matrix

i j r                
is close to the matrix

2 i (σ ) d iag with 2 2 i i a σ = 2
. Hence, the following approximation can be written:

2 i (σ ) T R d iag ≈ φ φ ≈ φ φ ≈ φ φ ≈ φ φ As 1 T M - -- - φ = φ φ = φ φ = φ φ = φ , (16) implies that 2 2 i i (σ ) (σ ) T 1 R M d ia g M d ia g - -- - ≈ φ φ = φ φ ≈ φ φ = φ φ ≈ φ φ = φ φ ≈ φ φ = φ φ
This shows that the columns of φ φ φ φ are close to the eigenvectors of the matrix R * M.

We will give sufficient conditions to ensure that the matrix

i j r                
is close to the matrix

2 i (σ ) d iag .
In the following sections, we first recall the classical condition for the sampling rate 1 ∆t , and then we study the sufficient conditions concerning the total observation time T m ∆ t = = = = for the case of an undamped system. Finally, we generalize these conditions for the case of a system with light proportional damping.

Choice of the sampling rate

We recall that the average of a continuous function g on an interval [0, T] is the limit for m → + ∞ of the average values taken at k t , with k t = k ∆ t and m ∆ t = T . This can be expressed as follows:

m k 0 k 1 1 1 g(t ) g(t) dt with m T t m t = T m t T → +∞ =    ∆ ∆   →  ∆   ∑ ∫ (22)
The necessary condition for

k 1 1 ( ) m k g t m =
∑ to be reasonably approximated by the average value of g on [0, T] is the verification of the Nyquist-Shannon sampling theorem, i.e., the sampling rate 1 ∆t should be higher than twice the maximum frequency of the signal g.

For the discrete dynamic system presented in (3), frequencies are given by

i i ω f = 2π
. Therefore, one must choose ∆t such that:

1 n 1 1 max( ) i i ω ∆t ≤ > π (23)
Then, for T = m ∆ t , where T is the observation time and m the number of samples, the approximation is as follows:

0 1 1 1 ( ) ( ) ( ) ( ) m T i j i k j k i j i j k r r t r t r r t r t dt m T = = ≈ = ∑ ∫ ɶ (24)

Choice of total observation time T

In this section, we study the limit of the matrix [ ] i j rɶ when T increases. We will first consider the undamped case (3.2.1), and then the case with proportional light damping (3.2.2).

Undamped system

Equation (13), used in the undamped case, gives the following expression of the temporal evolution of each mode:

i i ( ) c o s ( ) i i r t a t = ω -θ (25)
We will prove that when certain conditions are verified, the matrix [ ] i j rɶ tends to

2 i a 2       diag , when m → + ∞ , with: 0 1 ( ) ( ) dt. T i j i j r r t r t T = ∫ ɶ
In greater detail, for when i a and j a are nonzero, we will show that:

j 1 2 2 ij ii 2 i i r r and a a 2 a ε ε -≤ ≤ ɶ ɶ (26) with M T ε T =
, where the constant M T is function of the analyzed signal and T is the time of observation. By increasing T, ε might arbitrarily become small. In the following, we will consider the case where i = j , and then the case where i j ≠ .

• The diagonal terms are:

2 i i 0 1 cos ( ) T 2 i i i r a t dt T = ω -θ ∫ ɶ (27) As 2 i i i i i 0 i sin(2( )) sin(2( )) 1 1 cos ( ) 2 4 T T T t dt T - ω -θ + -θ ω -θ = + ω
∫ the following estimates hold: This shows that the first part of ( 26) is verified if max T T ≤ ε .

T 2 i i 0 i i 1 2 1 1 2 cos ( ) 2 4 2 4 t dt T T T - ≤ ω -θ ≤ + ω ω ∫ (28) Therefore, i i 1 1 1 1 2 2 2 min i i 2 i r T T a -≤ ≤ ω ω ɶ (
• The non-diagonal terms are:

i i j j 0 cos( ) cos( ) T i j i j a a r t t dt T = ω -θ ω -θ ∫ ɶ
For the sake of simplicity, we assume that the phase angles are zero. Therefore:

i j i j i j 0 i j i j 1 1 1 cos( ) cos( ) sin( ) sin( ) 2 T t t dt T T   ω ω = ω + ω + ω -ω   ω + ω ω -ω     ∫ which leads to: i j 0 i j i j 1 1 1 cos( ) cos( ) 2 T t t dt     ω ω ≤ +   ω + ω ω -ω   ∫ (31) Since
i j i 2 m in ( ) ω + ω ≥ ω , one has: i max i j i max( ) 1 1 T ( ) 2 min( ) 4 T
≤ = = ω + ω ω π (32) Defining gap
T by:

gap i j i j 1 1 T min 2 min f f = = ω -ω π - (33)
where i f are the modal frequencies, we obtain:

max gap T T 1 2 ij i j r a a T + ≤ ɶ
When two frequencies are close to each other, gap T can be quite large in respect to m a x T . Thus, in order to have a good approximation of the matrix

2 i a 2       diag
, one should make up for this by increasing the observation time T.

• To conclude, in the undamped case, the coefficients of the matrix [ ] i j rɶ are equal to those of the matrix

2 i a 2       diag with a relative error of order ε 2 for M T ε T = , where M m ax gap T T T = + (34)
Finally, let us note that if for a certain 0 i , 0 i a is null, 0 j i r ɶ and 0 0 i i r ɶ are also null, which means that the 0 i t h mode is not activated, thus it is not necessary to take it into account in the previous inequalities.

Proportionally damped system

We recall that for a system with proportional damping, ( ) i r t is given by equation (13). As for the undamped system, we simplify the presentation by assuming that all phase angles are null. Using the same procedure as in section (3.2.1), the aim is to show that, under certain conditions, the matrix [ ]

i j rɶ approaches 2 i a 2       diag
, when T is large enough. In greater detail, when i a and j a are nonzero, conditions are set in order to insure that

1 2 2 ii 2 i r a ε -≤ ɶ and ij i j r ε a a 2 ≤ ɶ
for ε arbitrarily small. The cases i = j and i j ≠ will be studied.

• The diagonal terms are:

i i 2 t 2 i 0 e cos ( ) 2 T i ii a r t dt T -ω ζ = ω ∫ ɶ ɶ (35)
As the Taylor series of x eis an alternating series for x 0 ≥ , one has we have:

i i 2 t i i i i 1 2 1 2 e 1 T t -ω ζ -ω ζ ≤ -ω ζ ≤ ≤ ɶ ɶ
Since the function given by

2 i c o s ( ) t ω ɶ
is positive, we deduce that:

i i 2 t 2 2 2 i i i i i 0 0 0 1 1 ( 1 2
) cos ( ) e cos ( ) cos ( )

T T T 1 T t dt t dt t dt T T T -ω ζ -ω ζ ω ≤ ω ≤ ω ∫ ∫ ∫ ɶ ɶ ɶ ɶ ɶ
Using the inequality (28) one has:

i i i i 1 1 1 1 (1 2 ) 2 4 2 4 ii 2 i r T T T a   - -ω ζ + ≤ ≤ +   ω ω   ɶ ɶ ɶ ɶ
which yields to:

i i i i i i i i i i i 1 1 1 (1 2 ) 4 2 4 T 1 1 1 1 4 2 2 4 ii 2 i ii 2 i r T T T a r T T T a -ω ζ --ω ζ ≤ - ≤ ω ω -ω ζ - -ζ ≤ - ≤ ω ω ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ
Using the following notations i max min max(T ) 1 T , 4 2 min( ) max( )

i i i 1 T and 0 = = = ≤ ζ ≤ η π ω ω ɶ ɶ ɶ ɶ ɶ one has max max i i i i i min T T 1 1 1 1 2 2 4 2 2 4 2 2T ii 2 i r T T T T T T a η -η - - ≤ -ω ζ - -ζ ≤ - ≤ ≤ ω ω ɶ ɶ ɶ ɶ ɶ ɶ ɶ (36)
If T is large enough and if the damping ratio η is low enough to satisfy the following inequalities:

max min T , T T ε ε an ε T d ≤ η ≤ η ≤ ɶ ɶ (37)
then the estimation (36) allows one to prove that the inequality max min T 3 T

T T η + + η ≤ ε ɶ
is verified, and so:

3 1 2 2 2 ii 2 i r ε ε a - ≤ - ≤ ɶ
• The non-diagonal terms are: i i j j ( ) t i i j j 0 e cos( ) cos( )

T i j i j a a r t t dt T -ω ζ -ω ζ = ω -θ ω -θ ∫ ɶ ɶ ɶ
The function

∫ ∫

As previously in (31):

c i i j j 0 i j i j 1 1 1 cos( ) cos( ) 2 t t dt     ω -θ ω -θ ≤ +   ω + ω ω -ω   ∫ then: 0 i j i j e 1 1 2 ij i j r a a T     ≤ +   ω + ω ω -ω   ɶ ɶ ɶ ɶ ɶ
which is the same result found in the undamped case. In order to have

2 ij i j r ε a a ≤ ɶ , with ε set in advance, it is sufficient to choose: M M max gap T T T T with T ≤ ε = + (38)
In conclusion, the conditions (37) and ( 38) are sufficient for the coefficients of the matrix [ ] i j rɶ to be approximations, of order 3 2 ε , of those of the matrix

2 i a 2       diag
. It is important to note that for a given damping ratio η , one cannot expect to have an approximation of degree ε smaller than max min T T η ɶ ɶ .

Numerical results and conclusion

To illustrate this theoretical study of the normal modes approximation accuracy by the POD, some numerical tests are presented in the following.

The relative error is calculated for different elementary initial data and for different observation times. The figures present the evolution of the mean relative errors with respect to T/T max .

To begin, we consider a diagonal matrix M. The example presented was previously chosen by Feeny and Kappagantu in [1 section(2.2)]. They used a total observation time T equal to four fundamental periods: . One can notice that the error is less than 10% when T is greater than m a x 1 0 T (which is about the same T considered in [1 section(2.2)], and that the accuracy is equal to 0.01589 and is coherent with Feeny's results. Figure [START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF] proves that, for systems with distinct modal frequencies ( To conclude, it is now numerically and theoretically proven that the POD, under the sufficient conditions discussed in this paper, can be used for dynamical low damping systems with non-diagonal M, even when two modal frequencies are close together. This careful analysis can be useful for finding the cause of poor approximation properties in more complex cases, such as analysis of nonlinear normal modes in vibrating systems.

  positive definite, their eigenvalues are positive and non null. Hence, one can define

Figure ( 1 )

 1 Figure[START_REF] Feeny | On The Physical Interpretation of Proper Orthogonal Modes in Vibrations[END_REF] illustrates the mean relative error evolution of the approximated normal modes, all the way to m a x 1 0 0 T T =

  is acceptable (lower than 10%) when we vary T/T M and do not consider gap T (since the gap between the frequencies is large).The next example concerns a non-diagonal, symmetric and positive definite mass matrix M, and two different matrices K, such that the modal pulsations of the system are respectively ω = [1, 2, 3] and ω = [1, 2, 2.(2) shows the poor accuracy obtained in the second case for the same T/T max . In the first case, the second one, since it has two similar modal frequencies to have a relative error of order epsilon, it is sufficient to have an observation time T less than * m ax T ε . But when two frequencies are too close to each other, in order to obtain an accuracy of order epsilon, it is necessary to take into account (3) considers the evolution of the relative error when T/T M increases for different damping ratios ζ . When M T T ≥ , one should examine the values of the relative error for three different values of proportional damping ζ . The comparison leads one to deduce that the proportional damping ζ must be low. Otherwise, the relative error becomes large and thus unacceptable.


A Dietz

C Julcour

A M Wilhelm

H Delmas

Selective hydrogenation in trickle-bed reactor. Experimental and modelling including partial wetting

Keywords: Trickle-bed reactor, selective hydrogenation, trickle-bed modelling

A steady state model of a trickle bed reactor is developed for the consecutive hydrogenation of 1,5,9-cyclododecatriene on a Pd/Al 2 O 3 catalyst. Various experiments have shown that the selectivity of this reaction towards the product of interest is much lower in co-current down-flow (trickle-bed) than in up-flow. This is due to uneven liquid distribution and to partial wetting of the catalyst surface at low liquid flow rates. The non-isothermal heterogeneous model proposed here takes into account the partial wetting of the catalyst, as well as the resistances to heat and mass transfer at the gas-liquid, liquid-solid and solid-gas interfaces. It assumes that the catalyst particles can be divided into two distinct concentration zones corresponding to the wetted and dry catalyst surfaces; mass transfer between these two zones is described by a simplified diffusion mechanism. Compared to previous models assuming a uniform concentration of liquid-phase components inside the catalyst particles, this model improves the prediction of the outlet concentrations of hydrogenation products.

Introduction

Fixed-bed reactors with down-flow of gas and liquid (trickle-bed reactors) have been extensively studied and widely used in industrial practice for many years. In recent works, the performance of such reactors was studied for various operating conditions and several models were proposed for the studied reactions.

Three phase catalytic reactions are most often carried out in trickle-bed reactors despite poor thermal control, uneven liquid distribution and, at low liquid flow rate, partial wetting of the catalyst. This phenomenon of partial wetting may have positive effects on the reaction rate when gas liquid mass transfer is the major rate limitation. Nevertheless it reduces selectivity as far as consecutive reactions are concerned. Models accounting for partial wetting are generally oversimplified due to two major problems. The first one arises through the poor knowledge of wetting, only quantified by an overall external wetting efficiency of the catalytic bed, which determination is highly questionable and which does not describe local features at the particle level : are all the pellets equally partially wetted or on the contrary does this wetting efficiency corresponds to the number fraction of pellets completely wetted, the remaining pellets being completely dry? The second point concerns the coupling of reaction and mass transfer inside the catalyst pores where complete internal wetting by capillarity is generally assumed. In the case of partially wetted pellets, the unsymmetrical conditions at the catalyst surface should lead to multidimensional diffusion equations which are rarely taken into account in the modelling. Simplified expressions of the effectiveness factor are usually used, accounting only for concentration gradients of the limiting reactant. [START_REF] Tan | Catalyst particle effectiveness with unsymmetrical boundary conditions[END_REF] compared those approximate solutions (weighting factor models, as reported in [START_REF] Ramachandran | Effectiveness factors in trickle-bed reactors[END_REF]) to rigorous results for slab geometry with two faces of equal areas covered at different reactant concentrations (wetting factor = 0.5) and cubic particle with one or more faces covered by liquid. They concluded that those simplified expressions give reasonably good results for a simple first-order (gas-limited) reaction. However concerning complex reactions and selectivity analysis only few works have been dedicated to the effects of partial wetting, although in any reactor modelling selectivity provides a very sensitive tool, with which no fitting parameter can be adjusted and may lead to possible misunderstanding of actual features. For wetting efficiency determinations a pseudo effective 1D diffusivity has sometimes been identified from tracer analysis or effective reaction rates but most often reactor models keep diffusion unchanged. In other words no comprehensive statement of the liquid reagent behaviour under the dry surface of the catalyst has been addressed up to now. Vergel (1993) studied the selective hydrogenation of butadiene in trickle-bed and flooded-bed reactors. He noticed that the selectivity remains stable in up-flow in a wide range of fluid velocities, but decreases in the trickle-bed at low velocities. The author proposed that the flow conditions inside the down-flow reactor can be represented by two parallel plug flow reactors of equal volume, receiving different fractions of the liquid and gas flow rates but completely wetted. [START_REF] Rajashekharam | A trickle-bed reactor model for hydrogenation of 2,4-dinitrotoluene: Experimental verification[END_REF] investigated the hydrogenation of 2,4-dinitroluene in a trickle-bed reactor. They combined the contributions of stagnant liquid covered as well as unwetted zones on the catalyst surface, but neglected liquid-solid and intraparticle mass-transfer resistances for the liquid-phase components. In this work the selective hydrogenation of 1,5,9-cyclododecatriene (CDT) to cyclododecene on a shell-type Pd/Al 2 O 3 catalyst was investigated for being a complex and highly exothermic reaction of industrial interest. The consecutive hydrogenation steps result in the formation of the following products: CDD (cyclododecadiene), CDE (cyclododecene), and CDA (cyclododecane). [START_REF] Julcour | Selective hydrogenation of 1,5,9-cyclododecatriene in up-and down-flow fixed-bed reactors: experimental observations and modeling[END_REF] worked on the same reaction and same reactor as in the present study (diluted bed). The authors observed a poor heat transfer and a low selectivity in the down-flow mode. A classical model assuming a uniform catalyst wetting but including direct gas-solid mass transfer to the dry catalyst fraction was proposed. The concentration gradients of the liquid-phase components inside the catalyst were also neglected. This model could fit the overall conversion, but not the selectivity due to heterogeneities of liquid flow and partial wetting. [START_REF] Chaudhari | Hydrogenation of 1,5,9-cyclododecatriene in fixed-bed reactors: down-vs. upflow modes[END_REF] investigated the same reaction in a smaller bench-scale reactor. In their short and non-diluted trickle-bed, channeling of liquid flow being not significant, the model of [START_REF] Rajashekharam | A trickle-bed reactor model for hydrogenation of 2,4-dinitrotoluene: Experimental verification[END_REF] satisfactorily described the behavior of the reactor.

Experimental set-up

The experimental set-up that was used is presented elsewhere: Julcour (1999) and [START_REF] Julcour | Selective hydrogenation of 1,5,9-cyclododecatriene in up-and down-flow fixed-bed reactors: experimental observations and modeling[END_REF]. The reactor consists of a jacketed packed-bed column with an inner diameter of 0.026 m and a bed height of 1.2 m. The tube is filled with a mixture of 1/6 of catalyst pellets (0.5% Pd on alumina from Degussa) and 5/6 of inert alumina pellets of the same shape. The liquid reactant is also diluted, as previous experiments have shown that hot spots could not be avoided otherwise (Julcour, 1999): a 15 wt % solution of CDT in isododecane is used as feedstock. Before entering the reactor gas and liquid phases are mixed in an inert fixed bed in order to homogenize the temperature and to achieve gas-liquid equilibrium.

Five automatic valves commanded simultaneously by a switch are placed along the gas and liquid circuits in order to manage the fluid circulation mode (up-and down-flow). Seven temperature probes and six liquid sampling valves are located along the reactor in order to measure axial temperature and concentration profiles. The control of the reactor temperature is achieved by a fast flowing of thermal oil (Marlotherm) in the jacket at a nearly constant temperature. Table 1. Operating conditions.

Modelling

Model (I): Uniform catalyst wetting, and uniform concentrations in liquid-phase components inside the catalyst (classical model)

First of all, the model used by [START_REF] Julcour | Selective hydrogenation of 1,5,9-cyclododecatriene in up-and down-flow fixed-bed reactors: experimental observations and modeling[END_REF] was applied. This model assumes that the surface of the catalyst pellet can be divided into two zones (fig. 1): a dry zone where the gaseous reactant is in direct contact with the catalyst, and a wetted zone where hydrogen is transferred through the liquid phase. The organic compounds are supplied to the catalyst though the wetted zone only.

Fig 1. Classical model of a partially wetted particle.

The assumptions of the model are the following ones:

-gas and liquid are in plug flow; -gas and liquid are considered as a pseudo homogeneous fluid for the energy balance; -the liquid-phase reactant is non-volatile; -internal diffusion limitations (intraparticle resistance for H 2 ) are taken into account through an apparent kinetics; -complete internal wetting is achieved due to capillary forces; -the catalyst temperature is uniform; -the contribution of stagnant liquid covered zones is negligible; -the liquid-phase reactant is in large excess with respect to the gaseous reactant so that concentration gradients of the liquid-phase components inside the catalyst may be ignored;

E Ex xt te er rn na al ll ly y d dr ry y z zo on ne e H 2 H 2 hydrocarbons W We et tt te ed d z zo on ne e

-severe intraparticle gradients are expected for H 2 which is consumed within a very thin zone close to the catalyst surface; -the catalyst wetted fraction is the same for all pellets (f w ).

The equations of the model are : Mass balances

• Liquid phase (plug flow) Hydrocarbons (k) ( ) ( ) k S, k L, k LS c w k L, C C k a f z F 0 - - ∂ ∂ - = (1) k = CDT,CDD, CDE, CDA H 2 ( ) ( ) ( ) H2 L, H2 L, L H2 Sw, H2 L, H2 LS c w H2 L, C C a k C C k a f z F 0 - + - - ∂ ∂ - = * (2) • Catalytic phase Hydrocarbons (k) ( ) ( ) ( ) ( ) ( ) ( ) S H2 Sd, m S, d i, w S H2 Sw, m S, w i, w i ik B k S, k L, k LS c w T , C , C ' r f 1 T , C , C ' r f C C k a f 0 - + ν ρ + - = ∑ (3) k = CDT, CDD, CDE, CDA m = {CDT, CDD, CDE} H 2 -wetted part ( ) ( ) ( ) ∑ ν ρ + - = i w i, w iH2 B H2 Sw, H2 L, H2 LS c w ' r f C C k a f 0 (4) -dry part ( ) ( ) ( ) ( ) ( ) d i, w i iH2 B H2 Sd, H2 L, H2 GS c w ' r f 1 C C k a f - 1 0 - ν ρ + - = ∑ * (5) Energy balances • Fluid phase (gas + liquid) ( ) ( ) ( ) ( ) W1 F R 1 w F S FS c F T T V A h T T h a z Fh 0 - - - + ∂ ∂ - = (6)
• Catalyst phase

( ) ( ) ( ) d i, w w i, w i Ri B F S FS c ' r f 1 ' r f H T T h a 0 - + ∆ ρ - - - = ∑ (7) • Cooling oil (Marlotherm) ( ) ( ) ( ) M W1 an 1 2 M T T V A h Fh z 0 - + ∂ ∂ - = ( ( ) ( ) M W1 1 2 W1 F 1 w T T A h T T A h 0 - - - = ) (8)
5/12

Model parameters

Model (I) involves around 20 parameters representing complex kinetics, heat and mass transfer.

The parameters of the apparent kinetics were identified form batch experiments carried out with the 3.1 mm diameter pellets.

A Eley-Rideal model was selected:

∑ = + = α 3 1 j j j H j j i i C K 1 C C K k ' r i 2
Parameters can be found in Julcour (1999).

(

) 9
where i is the reaction number, and j is the component number ( j = 1: CDT, 2: CDD, 3: CDE). This model has been first assessed for up-flow mode, using a wetting efficiency equal to unity and empirical correlations established by Stüber (1995). It can conveniently predict the axial concentration and temperature profiles (cf. figure 2). (where the relative molar fraction (or %) of CDT is given as the ratio of CDT concentration to the sum of CDT, CDD, CDE and CDA concentrations).

In down-flow mode, coefficients k L a, k LS for H 2 , k LS for hydrocarbons, k GS , h FS , and h w are calculated by correlations from the literature [START_REF] Fukushima | Liquid phase volumetric and mass-transfer coefficient, and boundary of hydrodynamic flow region in packed column with current downward flow[END_REF], [START_REF] Tan | A dynamic method for liquid-particle mass transfer in trickle beds[END_REF], [START_REF] Satterfield | Liquid-solid mass transfer in packed beds with downward cocurrent gas-liquid flow[END_REF], [START_REF] Dwidedi | Particle-fluid mass transfer in fixed and fluidized beds[END_REF], [START_REF] Whitaker | Forced convection heat transfer correlations for fluid in pipes, past flat plates, single cylinders, simple spheres, and for fluid in packed beds and tube bundles[END_REF], [START_REF] Specchia | Heat transfer in trickle-bed reactors[END_REF]). Only the wetting efficiency (f w ) has been used as a fitting parameter. In the investigated range of fluid velocities, the correlation of [START_REF] El-Hisnawi | Trickle-bed reactors: Dynamic tracer tests, reaction studies, and modeling of reactor performance[END_REF] gives a value around 0.7 for the catalyst wetting efficiency, but as shown in fig. 3, a small deviation in the estimation of this factor results in dramatically different predictions of the outlet concentrations. However, in the particular case investigated fig. 3, this parameter f w was found to have a low influence on the concentration of CDE. In fact as severe hydrogen limitations are governing the reaction rates increasing the dry fraction accelerates similarly all reactions. In any case CDT concentration should decrease and CDA increase, but the behaviour of intermediates depends on the overall hydrogenation advancement. As the same catalytic bed was used for many runs, a deactivation factor (deact) is also introduced in the model to correct the rate equations measured with fresh catalyst (r' i,actual = deact × r' i ). However the catalyst activity was found to decrease very steeply during the first runs after its loading and then only very gradually (Stüber, 1995). The four previous (consecutive) experiments were performed in the latter conditions, so that a constant deactivation factor could be used: deact = 0.85.

Fig 3 . Influence of the catalyst wetting efficiency.

(where reaction advancement is calculated by: 3

3.X 2.X X CDA CDE CDD H2 + + = Ω
, while selectivity S is defined as the ratio the considered product to the sum of all the reaction products (CDD, CDE, CDA)).

Results

Based on the following criteria: The catalyst wetting efficiency is approximately 0.7, as expected by the correlation of El-Hisnawi. Nevertheless, in the case of high gas flow rate, predictions can meet experimental hydrogenation rates only for a lower value. This can be explained by the fact that the correlation for f w does not account for gas flow rate variations.

( ) 2 ∑ - = i iC iE X X C (i =
Even though the model can predict the overall hydrogenation rate, it can never fit the selectivity values for CDE (see fig. 4). One reason may be that the model assumes the hydrocarbon concentrations inside the catalyst pellets to be uniform. In fact, it seems obvious that the concentrations of the liquid-phase components are not the same everywhere inside the catalyst because hydrogen transfer is greater to the dry surface than to the wetted surface and a diluted reactant is used. This can have in return some effects on the reaction rates, especially when a zero order with respect to hydrocarbons is not still verified.

Selectivity (S)

S(CDD) S(CDE) S(CDA) Ω H2

It is then necessary to differentiate inside the catalyst particle two concentration zones for hydrocarbons corresponding respectively to the wetted and dry surfaces. Moreover in model (I) the wetting efficiency is supposed to be the same for all the pellets, while in fact channeling of liquid flow has every chance to occur in this long reactor at low liquid velocities, leading to an heterogeneity of the wetting at the bed scale.

Model (II): Reactor split into two portions operating in parallel under different flow conditions

To account for wetting heterogeneitites at the reactor scale, a very simple model was derived from model (I), assuming the reactor to be a linear combination of two reactors working in parallel under different wetting conditions: one (volume fraction k) in which the catalyst is totally wetted (wetting efficiency = 1) and the other one (volume fraction (1-k)) in which the wetting efficiency value is ranging between 0.5 and 0.8. After a trial and error procedure, the following results have been obtained, showing that model (II) predictions can be closer to experimental data, especially when regarding the outlet concentration of the final product CDA. This is due to lower values of the wetting efficiency used in the partially wetted section: in this part, the reaction proceeds quickly leading to complete hydrogenation. However the gaps for CDD and CDE are still significant (fig. 5).

Model (III): Uniform catalyst wetting, concentration heterogeneities in the particle for liquid-phase compounds

As the previous models could not predict the selectivity values, a new model was derived to improve the description of the catalyst scale phenomena under partial wetting conditions. As the first model, model (III) assumes that the external wetted fraction is the same for all the pellets, but it accounts for concentration heterogeneities of the liquid-phase components inside the catalyst: each particle is divided into two distinct concentration zones corresponding respectively to the wetted and dry catalyst surfaces (see fig. 6).

Hydrocarbons are supplied to the catalyst through the wetted surface only, while hydrogen passes both through the liquid phase and directly at the dry surface. Mass transfer between these two zones is described by a simplified diffusion mechanism, using arbitrarily the dry surface as a reference:

( ) ( )

i Sd i Sw cd Li c w i C C L D a f F , , 1 - - = for i = CDT, CDD, CDE, CDA, H 2
where L cd is the diffusional characteristic length. These considerations lead to modify equations ( 3), ( 4) and ( 5): two zones (wetted and dry) and thus two sets of concentration are now considered for hydrocarbons and the exchange term between the two zones is also added.

Influence of model parameters

This model involves one more parameter than model (I): the diffusion characteristic length L cd .

As in model (I), the wetting factor has a significant effect on the overall hydrogenation rate and on the selectivity. The increase of the diffusional characteristic length reduces the reaction rate because it reduces the hydrocarbon fluxes towards the catalyst dry zone. This parameter has thus a major effect on selectivity (fig. 7).

H 2 H 2 hydrocarbons

9/12 Fig. 7. Influence of the diffusional characteristic length on the product distribution. As expected, the wetting efficiency was found to increase with the liquid flow rate, and less obviously to decrease with the gas flow rate. Its order of magnitude is near to the one predicted by the correlation of El-Hisnawi. For a constant liquid velocity, the diffusional characteristic length decreases when the gas flow rate increases. An increase of the gas flow rate should result in a better spreading of the liquid at the catalyst surface so that the catalyst dry fraction would not correspond to a single region, but to several ones. The real exchange surface between dry and wetted zones increasing, this leads in the model to a decrease of the diffusional characteristic length.

Optimal parameters

On the other hand increasing the liquid flow rate would result in less dry zones then in higher diffusional characteristic length.

Results

This model predicts the molar fractions of CDE and CDA correctly, but there are still small differences for CDT and CDD. Nevertheless, this model predictions are very close to the experimental values (fig. 8).

Fig 2 .

 2 Fig 2. Axial concentration profiles in up-flow mode. (P = 4 bars, T Wall = 100°C, u Le = 0.53×10 -3 m/s, u Ge = 0.042 m/s).

Fig 4 .

 4 Fig 4. Distribution of products obtained experimentally and predicted by model (I).

Fig 5 .

 5 Fig 5. Results of the model with two parallel reactors.

Fig 6 .

 6 Fig 6. New model involving different hydrocarbon concentrations in the two zones (wetted and dry).

  (based on the sum of the four hydrocarbons: CDT, CDD, CDE, CDA) z = reactor length, m ∆H R = heat of reaction, J/mol ρ B = density of catalyst bed, kg/m 3 ν = stoichiometric coefficient Ω H2 = overall hydrogenation advancement Subscripts and abbreviations: C = calculated; d = dry; E = experimental; F = fluid; G = gas; L = liquid; M = Marlotherm; S = solid; W = wetted; W1= inner reactor wall.

Table 2 .

 2 Identified values for wetting efficiency.

  Table3shows the optimised values of the wetting efficiency and diffusional characteristic length for the new model.

		fw	Lcd (mm)	Criteria
	Exp 1	0,65	1,16	0,0078
	Exp 2	0,55	0,44	0,0034
	Exp 3	0,6	0,96	0,0044
	Exp 4	0,65	1,20	0,0038

Table 3 .

 3 Optimised values for wetting efficiency and diffusional characteristic length.

Comparison of the different models

The following graphs (fig. 9) compare the different models for each experiment. The model with two parallel reactors always predicts the fraction of CDT and CDA at the outlet of the reactor correctly, but model (III) fits the experimental results more precisely.