M Turki

M Kebdani

G Dauphin-Tanguy
email: genevieve.dauphin-tanguy@ec-lille.fr.

A Dazin

P Dupont

Experimentally Validated Bond Graph Model of a Brazed-Plate Heat Exchanger (BPHE)

Keywords: Fluid loop, heat exchanger, bond graph, dynamic, single phase flow

The paper deals with the Bond Graph (BG) modeling and the model validation of a brazed-plate heat exchanger. This device is an important part of a mechanically pumped cooling loop. A thermo hydraulic BG model is developed and compared with experimental data. Optimization is performed to determine the best value of the convection heat exchange coefficients to be fixed in the model.

I.

INTRODUCTION

Power components dissipate heat flows which represent a significant source of heat. If it is not processed properly, this results in problems that lead to a complete dysfunction of the component. Several solutions have been proposed to prevent the damage of embedded devices and ensure their proper functioning.

Based on these solutions, the choice of the enterprises has been oriented towards mechanically pumped cooling loops that are undoubtedly the most efficient in terms of heat transfer (Kebdani et al [START_REF] Kebdani | Dupont Bond Graph Model of a mechanically Pumped Biphasic Loop[END_REF]). Indeed, in addition to their architectural flexibility, they guarantee appropriate cooling in suitable operational conditions. It is a very important property especially when we know that the continuous cycling of the temperature has a negative impact on the electronic components' life.

The figure 1 shows up the considered cooling loop. It is composed of a pump, a pressure regulator, an evaporator, a condenser and pipes.

Various heat exchanger technologies exist; the choice depends on the nature of intended use. For example, for space activities, one can use a radiator; for land-based activities (automotive, rail...) an air exchanger with air crossflow may be adequate.

Fig.1: Design of the cooling loop

A single-phase fluid loop is a heat transfer circuit. It is generally a closed loop, wherein a fluid is initially in the liquid state. The fluid is heated to store thermal flux from a heat source, and is then transferred to a heat sink where the hot fluid is cooled and even sub cooled releasing heat to a cold source. Such fluidic loops have high cooling efficiency.

For our purpose, we have chosen to work with the condenser SWEP of reference: "B5Tx6" (Fig. 2

Fig.2: Design of the condenser

Also, the corrugated structure of the thermal plates and chevrons promote turbulent flow (Focke et al. [5]). This leads to high heat transfer efficiency and consequently better performance of the whole installation. This paper is organized into five sections; the first section is a presentation of our work's overall context and focuses on what has already been done in the BG condenser modeling. The second part deals with the BG model proposed by introducing various assumptions taken into account and the equations put in. The third section is dedicated to the experimental set-up and the various tests performed. The fourth section illustrates an optimization of the heat transfer coefficient. The last part concerns the validation of the proposed model with experimental results

State of the art

The evaluation of such exchanger's performance is not easy as it needs the development of specific experimental methods, especially for prediction of the exchange coefficient and pressure drops. In this context, various authors conducted studies of instrumentation and visualization:  Among the works related to this type of heat exchangers, there have been attempts based on the observation of twophase water-air flows (Vlasogiannis et al.  Xiaoyang et al [START_REF] Rong | Two-phase header flow distribution in a stacked plate heat exchanger[END_REF] declare that it is possible to estimate theoretically the performance of a BPHE operating in single phase with water as refrigerant fluid. The final results of their work show that it would be appropriate to utilize correlations of pressure drop and heat transfer that take into account corrugation chevron angles.

 Performance of plate heat exchangers under single phase operations are extensively inspected since the 1980s (Cooper and Usher, 1983 [10]; Raju and Bansal, 1983 [11]; Focke et al., 1985 [12]; Shah and Focke, 1988 [13]; Bansal and Muller-Steinhagen, 1993 [14]). According to these authors, these theoretical works provide a reliable basis that ensures an efficient sizing of the exchanger.

 B. Ould Bouamama (1997) [START_REF] Ould Bouamama | Bond Graph modelisation of steam condenser[END_REF] proposes a model, based on BG methodology, of a simple tubular condenser. The mathematical formulation of the problem is clearly described and used to generate a BG intending to predict the dynamic behavior of the condenser operating under two-phase conditions.

In our knowledge, there are no other attempts in the published literature which aim to model BPHEs by referring to the BG approach.

The main purpose of this paper is to propose a validated thermo hydraulic BG model of a chevron type BPHE, operating under single-phase conditions, where the dynamic regime is considered. Simulation results are validated using the test rig developed by the French enterprise Atmostat.

Moreover, according to the research mentioned above, the correlations of the heat exchange are strongly related to experimental conditions in which they were developed. This means that for our case it would be legitimate to pick up the most suitable correlation and, thereafter, perform an optimization study in order to minimize the difference between experimental results and those provided by the model.

Given the nature of the device which involves several physical domains (thermal, hydraulic, electric ...), the bond graph methodology appeared as an appropriate modeling tool.

II. Bond graph model of the condenser

The thermal fluid processes involve two phenomena: convection, describing the heat transfer between the fluid in motion and walls, and heat conduction, representing the thermal power transmitted due to a temperature difference. Furthermore, the friction effect in the fluid in movement contributes to the heating of the fluid, especially in case of complicated geometries. Then, it would be wise to analyze the thermo hydraulic coupling for a better modeling of the brazed-plate heat exchanger. This is the purpose of the present section.

1) Generalized variables used.

 Hydraulic part: Pressure and mass flow rate for which the product is not a power. The developed model is a "pseudo bond graph".

 Thermal part: for a fluid in motion the entropy balance is not conservative, which justifies the necessity to choose the heat balance instead.

2) Coupled Bond Graph of the condenser.

The condenser is mainly composed of a (hot) primary circuit (PC) and a (cold) secondary circuit (SC) separated by brazed-plates as shown Fig. 3.

Assumptions:

1) The internal geometry of the brazed-plate condenser is too complicated to be precisely modeled. Thus we simply consider that such exchanger has a simple rectangular section ) * ( e e without any corrugations.

Fig.3: longitudinal cup of the condenser

2) The upstream pressure variation is known and modeled as MSe: effort source Fig. 4.

3) The device is not correctly isolated from the ambient, which is modeled as a MSe. 4) The flow is supposed to be constantly monophasic. 5) Working fluid and secondary fluid are pure water. 6) Each part of the condenser (PC and SC) is modeled as a unique volume where the phenomena are supposed to be homogeneous.

The figure 4 shows up the bond graph model of the condenser.

I. RC-elements

The hydraulic power transfer is represented Fig. 4

   (1)
where the Darcy coefficient « f » depends on the flow regime as shown in table.1.

] * ² * * 2 1 [ * 2 2   v P   (2)
For an elbow with sharp angle "90°" one can take ζ= 1.3.

Fig.4 Pseudo Bond Graph model of the condenser.

Thermal part:

 The heat flow due to the total pressure losses lead to

 m P P Q   * ) ( 2 1    

II. R-elements

The two elements R: conv deal with the heat exchange by convection between each working fluid and the brazed-plate that separate the two fluids. The corresponding heat flow is:

T S h Q conv      (6)
The two elements R: cond_ext represent the heat transfer by conduction within the wall of the condenser to the ambient. The corresponding heat flow is:

T S K Q equiv cond      (7)
The problem is now to experimentally determine the values of the convection coefficients PC h and SC h .

CS CP cond equiv

h h h K 1 1 1 1    (8)

III. Experimental set-up

The set-up presented in Fig. 5 has been designed by the French company Atmostat. It is composed of all the devices of the cooling loop. The test bench was equipped with pressure, temperature, and flow rate sensors. All the experiments have been realized in single-phase state of the fluids. In fact, it can be argued that as long as the flow is purely monophasic, the heating Q  due to the friction may be omitted as shown in Fig. 7 where wall temperatures are exactly the same with and without thermo hydraulic coupling in PC, Fig. 7 (a) and in SC, Fig. 7 (b).

However, a biphasic study of the condenser will need to take into account the thermal-hydraulic coupling.

Comparison between fluid and wall temperature:

At the steady-state the difference between the core working fluid temperature (model), Fig. 8 (a

IV. Identification and optimization of the convective conductance.

The accurate determination of the heat transfer coefficient "h" inside the condenser is very difficult because of the complexity of the corrugated geometry of the exchanger. However the dynamic BG model proposed in this study is equipped with a semi empirical model for its estimation. The objective of the current section is to optimize the thermal conductance in both compartments of the condenser PC and SC. Starting with the nominal value " 0 h "calculated by the proposed semi empirical model (eq.9) for each compartment, then an optimization is performed on " 0 h " to minimize the difference in temperatures between the model and experience.

The optimization used in this study is based on Broydon Fletcher Goldfarb Shanno method, which is already integrated into the 20sim simulation software. This method uses both the gradient of a function and the second order gradient to determine the search direction. The search direction is kept for each new step until a minimum has been found. Then a new search direction is determined and the process goes on. a) First test with a constant conductance value. The brazed-plate heat exchanger used in our study is delivered with data sheet where the manufacturer mentions a value of h = 6530 C w/m² . The idea here is to launch the simulation of the first test case (4.a) with this imposed value and compare the profile of simulating temperatures with experimental measurements (Fig. 9).

The dynamic model simulated with this value of "h" shows a discrepancy of 2°C with the measured temperature in steady state. This means that the imposed value of conductance does not correspond to the right value. a) Improvement of the convective conductance. The scientific literature describes very few models of the convective conductance specific to brazed-plate heat exchanger. However, Alfa Laval is a condenser manufacturer which discloses a model adapted to our current application and whose formula is: 10 shows a better concordance between numerical resolution (red and blue curves) and experience (purple curves). The difference being recorded is reduced from 2°C to 1.2 °C in primary circuit, and from 0.9°C to 0.2°C in secondary circuit, which is inside the precision domain of the temperature sensors. The result is slightly better for the SC because the fluid is pure water, even though the fluid in the hot circuit (PC) is not really pure water but contains a certain quality of oil for the pump lubrication.

3274 , 0 ) ( 2 3 / 1    P P K h r       (

Test 2:

Initial conditions:  Cold mass flow rate = 4 cm3/s.  Hot mass flow rate = 0 cm3/s. At the hot circuit, we notice that at the steady state regime, the temperature difference between the curve from the model and the experimental is around 1 degree, which may correspond to the uncertainty of the thermocouples. Whereas

  ). See Condensers technical specification [2]. The BPHEs have been used for the first time in the 1930s. They were mostly integrated as mono-phasic (liquid-toliquid) exchangers in the food industry [3] and in the heat pumps [4] thanks to their multiple benefits listed below:  Optimized effective exchange surface.  Rational distribution of the flow in the channels.  Good adaptability.  High compacity.

[ 6 ]

 6 and Volker & Kabelac [7]) for a plate heat exchanger formed of a single channel constituted with a transparent plate. Consequently a flow pattern map is performed.  More recently Freund & Kabelac [8] have developed an experimental technique based on infrared visualization, to characterize the spatial distribution of the convective heat exchange coefficient for a single-phase flow (water).

  by blue half arrows, while the thermal part appears in orange half arrows. The coupling element RC represents both the pressure losses generated by walls friction (R part) and the thermal energy storage phenomenon in the exchanger (thermal capacitance effect). The equations describing these phenomena are: Hydraulic part: For this part, two kinds of pressure losses are taken into account in the dynamic model:  Linear pressure drop:

Table 1 :

 1 Darcy coefficient as a function of the flow regime. Singular pressure losses due to two elbows:



  Temperature in each part of the condenser is given by

Fig. 5 : 6 :Fig. 6 (

 566 Fig.5: Design of the real loop a) Hydraulic analysis.In view of validating the thermo-hydraulic model of the condenser, we consider here the following monophasic test:

  Fig.6 (a) Time evolution of the outlet pressure of the condenser. (b) Discrepancy between experimental pressure drop and modeled one.

Fig. 7 :

 7 Fig.7: Temperatures; in the PC (Graph a), in the SC (Graph b), time evolution of input flow in PC and SC (Graph c).

Fig. 8 :

 8 Fig.8: Time evolution of fluid temperature in PC (a), in SC (b), and flow evolution in both PC and SC.

Fig. 9 :

 9 Fig.9: Temperatures profiles in both PC and SC with SC PC h h  = 6530 C w/m² .

  K.s} Hdot_e_SC {J/K.s} [5] Focke, W.W., J. Zachariades, and I. Olivier. 1985. Effect of the corrugation inclination angle on the thermohydraulic performance of plate heat exchangers. International Journal of Heat and Mass Transfer 28:1469-1479.

  K.s} Hdot_e_SC {J/K.s}

Acknowledgment

This paper describes results from research supported by the FUI 14 Program in the context of project ThermoFluid RT labeled by the competitivity pole ASTech. The authors gratefully acknowledge Mr R. Albach from Atmostat company for his important contribution to the experimentation phase, and the members of the consortium for the multiple fruitful scientific discussions about heat exchanges.

at cold circuit, we clearly see that the two curves are superimposed (fig. 12).

Fig.10: Temperatures profiles in both PC and SC with

Optimized Alfa Laval "h".

V-Conclusion

In this paper we propose a Brazed-Plate heat exchanger model based on BG approach, taking into consideration the hydro-thermal exchanges that occur in the system. Due to simplified assumptions, the difference between theoretical and experimental results can be attributed to the uncertainty of the sensors. The model can be considered as reliable enough to represent the heat transfer in the heat exchanger in monophasic behavior. The model is parametrized in terms of device geometry, type of fluid, input variables (hot and cold sources and temperatures). The future research work is to model the two-phase regime and validate it through experiments.

Model variables.


Evgeny R Gafarov

Alexandre Dolgui
email: dolgui@emse.fr

Alexander A Lazarev

Frank Werner
email: frank.werner@ovgu.de

A New Effective Dynamic Program for an Investment Optimization Problem

Keywords: Multi-choice knapsack problem, Investment problem, Dynamic programming, FPTAS MSC: 90B30, 90C39, 90-08

After a series of publications of T.E. O'Neil et al. (e.g. in 2010), dynamic programming seems to be the most promising way to solve knapsack problems. Some techniques are known to make dynamic programming algorithms (DPA) faster. One of them is the graphical method that deals with piecewise linear Bellman functions. For some problems, it was previously shown that the graphical algorithm has a smaller running time in comparison with the classical DPA and also some other advantages. In this paper, an exact graphical algorithm (GrA) and a fully polynomial-time approximation scheme based on it are presented for an investment optimization problem having the best known running time. The algorithms are based on new Bellman functional equations and a new way of implementing the GrA.

Introduction

The selection of projects is an important problem in science and engineering. In general, project portfolio selection problems are hard combinatorial problems, often formulated as knapsack or in a more general form as multi-dimensional multi-choice knapsack problems, possibly also with multiple objectives. They have a wide range of business applications in capital budgeting and production planning.

Different techniques have been used to solve the knapsack problem or its generalizations. Since the project data are often uncertain or imprecise, sometimes fuzzy techniques are used. For instance, in [START_REF] Tavana | A fuzzy multidimensional multiple-choice model for project portfolio selection using an evolutionary algorithm[END_REF], a fuzzy multi-dimensional multiple-choice knapsack problem is formulated for the project selection problem, and then an efficient epsilon-constraint method and an multi-objective evolutionary algorithm are applied. In [START_REF] Chang | A fuzzy DEA and knapsack formulation integrated model for project selection[END_REF], a data envelope analysis, knapsack formulation and fuzzy set theory integrated model was suggested. In [START_REF] Beaujon | Balancing and optimizing a portfolio of R & D projects[END_REF], the problem of selecting projects to be included in an R&D portfolio has also been formulated as a multi-dimensional knapsack problem. If partial funding and implementation is allowed, linear programming can be used and the sensitivity of the project selection decisions is examined. The selection among ranked projects under segmentation, policy and logical constraints was discussed in [START_REF] Mavrotas | Selection among ranked projects under segmentation, policy and logical constraints[END_REF]. After ranking the projects by a multi-criteria approach, integer programming was applied to get a final solution satisfying the constraints. At the integer programming phase, a knapsack formulation was applied. Dynamic order acceptance and capacity planning on a single bottleneck resource has been considered in [START_REF] Herbots | Dynamic order acceptance and capacity planning on a single bottleneck resource[END_REF]. Stochastic dynamic programming was applied to determine a profit threshold for the accept/reject decision and to allocate a single bottleneck resource to the accepted projects with the objective to maximize expected revenue.

Often dynamic programming algorithms (DPA) are used for the exact solution and fully-polynomial time approximation schemes (FPTAS) based on these algorithms are derived for the approximate solution of such problems. In [START_REF] Kellerer | Improved dynamic programming in connection with a FPTAS for the knapsack problem[END_REF], a vector merging problem was considered in a dynamic programming context which can be incorporated into an FPTAS for the knapsack problem. Approximation algorithms for knapsack problems with cardinality constraints, where an upper bound is imposed on the number of items that can be selected, were given in [START_REF] Caprara | Approximation algorithms for knapsack problems with cardinality constraints[END_REF]. Improved algorithms for this problem were presented in [START_REF] Mastrolilli | Hybrid rounding techniques for knapsack problems[END_REF], where hybrid rounding techniques were applied. An efficient FPTAS for the multiobjective knapsack problem was given in [START_REF] Bazgan | Implementing an efficient FPTAS for the 0-1 multi-objective knapsack problem[END_REF]. It uses general techniques such as e.g. dominance relations in dynamic programming. Approximation algorithms for knapsack problems with sigmoid utilities were derived in [START_REF] Sristava | Knapsack problems with sigmoid utilities: Approximation algorithms via hybrid optimization[END_REF]. The authors combined algorithms from discrete optimization with algorithms from continuous optimization. In [START_REF] Guler | Algorithms with guarantee value for knapsack problems[END_REF], greedy algorithms for the knapsack problem were considered and improved approximation ratios for different variants of the problem were derived. Some applications of AI techniques to generation planning and investment were described in [START_REF] Wu | Applications of AI techniques to generation planning and investment[END_REF].

This paper deals with a project selection (or allocation) problem with the single criterion of maximizing the total profit under one budget constraint, which is also a generalization of the knapsack problem. More precisely, a set N = {1, 2, . . . , n} of n potential projects and an investment budget (amount) A > 0, A ∈ Z, are given. For each project j ∈ N , a profit function

f j (x), x ∈ [0, A],
is also given, where the value f j (x ) denotes the profit received if the amount

x is invested into the project j. The objective is to determine an amount x j ∈ [0, A], x j ∈ Z, for each project j ∈ N such that n j:=1 x j ≤ A and the total profit n j:=1 f j (x j ) is maximized. Closely related problems do not only exist in the area of project investment [START_REF] Li | Expanded Model of the Project Portfolio Selection Problem with Divisibility, Time Profile Factors and Cardinality Constraints[END_REF], but also in warehousing, economic lot sizing, etc. [START_REF] Dolgui | Supply Chain Engineering: Useful Methods and Techniques[END_REF][START_REF] Schemeleva | Lot-sizing on a Single Imperfect Machine: ILP Models and FPTAS Extensions[END_REF].

In the following, we work with piecewise linear functions f j (x). The interval [0, A] can be written as a union of intervals in the form

[0, A] = [t 0 j , t 1 j ] (t 1 j , t 2 j ] . . . (t k-1 j , t k j ] . . . (t kj -1 j , t kj j ]
such that the profit function has the form

f j (x) = b k j + u k j (x -t k-1 j ) for x ∈ (t k-1 j , t k j ],
where k is the number of the interval, b k j is the value of the function at the beginning of the interval, and u k j is the slope of the function. Without loss of generality, assume that b

1 j ≤ b 2 j ≤ . . . ≤ b kj j , t k j ∈ Z, j ∈ N, k = 1, 2, .
. . , k j , and that t kj j = A, j = 1, 2, . . . , n.

A special case of the problem under consideration is similar to the well-known bounded knapsack problem: maximize n j:=1 p j x j s.t.

n j:=1 w j x j ≤ A, x j ∈ [0, b j ], x j ∈ Z, j = 1, 2, . . . , n, (1)
for which a dynamic programming algorithm (DPA) of time complexity O(nA) is known [START_REF] Kellerer | Knapsack Problems[END_REF]. Dynamic programming algorithms and an FPTAS for the bounded set-up knapsack problem, which is a generalization of the bounded knapsack problem in which each item type has a set-up weight and a set-up value included into the objective function, were suggested in [START_REF] Mclay | Algorithms for the bounded set-up knapsack problem[END_REF]. These algorithms can also be applied to the bounded knapsack problem, and one of the dynamic programming algorithms has the same time complexity as the best known algorithm for the bounded knapsack problem.

It is known [START_REF] Posypkin | Speedup Estimates for Some Variants of the Parallel Implementations of the Branch-and-Bound Method[END_REF] that all branch and bound (B&B) algorithms with a lower and an upper bound calculated in polynomial time have an exponential running time close to 2 O(x) operations, unlike P = N P , where x is the input length. For example, for the one-dimensional knapsack problem, the time complexity is equal to or greater than 3 2

2 n+3/2 √ π(n+1)
, where n is the number of items. This means that B&B algorithms are not far away from a complete enumeration.

In a series of publications by T.E. O'Neil et al. [START_REF] O'neil | A Simple 2 O( √ x) Algorithm for PARTITION and SUBSET SUM[END_REF], it was shown that some of the well-known NP-hard problems can be solved by dynamic programming in sub-exponential time, i.e., in 2 O( √ x) operations. So, at the moment, dynamic programming seems to be the most promising way to solve knapsack problems. Some techniques are known to make dynamic programming faster, e.g., in [START_REF] Bar-Noy | Online Dynamic Programming Speedups[END_REF][START_REF] Eppstein | Speeding up Dynamic Programming[END_REF][START_REF] Wagelmans | Improved Dynamic Programs for Some Batching Problems Involving the Maximum Lateness Criterion[END_REF], functional equations and techniques were considered that are different from the ones considered in this paper.

The following problem is also similar to the problem under consideration:

minimize n j:=1 f j (x j ) s.t. n j:=1 x j ≥ A, x j ∈ [0, A], x j ∈ Z, j = 1, 2, . . . , n, (2)
where f j (x j ) are piecewise linear as well. For this problem, a DPA with a running time of O( k j A) [START_REF] Shaw | An Algorithm for Single-Item Capacitated Economic Lot Sizing with Piecewise Linear Production Costs and General Holding Costs[END_REF] and an FPTAS with a running time of O(( k j ) 3 /ε) [START_REF] Kameshwaran | Nonconvex Piecewise Linear Knapsack Problems[END_REF] are known.

In this paper, we present an alternative DPA based on a so-called graphical approach. This algorithm has a running time of O( k j min{A, F * }), where F * is the optimal objective function value. Thus, it outperforms an algorithm from [START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF] which has a worse running time close to O(nk max Alog(k max A)), where k max = max j-1,...,n k j . The second contribution of this paper is an FPTAS derived by a scaling argument from this new DPA. Note that an FPTAS was already proposed for the treated problem in [START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF], but the new FPTAS has an improved running time of O( k j n log log n/ε).

While the running time of a similar DPA from [START_REF] Shaw | An Algorithm for Single-Item Capacitated Economic Lot Sizing with Piecewise Linear Production Costs and General Holding Costs[END_REF] is proportional to the sum of all linear profit pieces times the budget A, the new algorithm replaces A by the largest profit of a single project. The main idea is as follows. Instead of evaluating the dynamic programming functions for every budget value t = 1, . . . , A, we keep the profit functions (depending on the budget value) as piecewise linear functions. These functions can be represented by a collection of linear pieces, instead of a full table of the values. Since all relevant data are integer and the profit functions must be non-decreasing in the budget value, one can easily bound the number of relevant linear pieces to obtain the improved complexity.

The remainder of the paper is as follows. In Section 2, we present the Bellman equations to solve the problem under consideration. In Section 3, an exact graphical algorithm (GrA) based on an idea from [START_REF] Lazarev | A Graphical Realization of the Dynamic Programming Method for Solving NP-hard Problems[END_REF] is presented. In Section 4, an FPTAS based on this GrA is derived.

Dynamic programming algorithm

In this section, we present a DPA for the problem considered. For any project j and any state t ∈ [0, A], we define F j (t) as the maximal profit incurred for the projects 1, 2, . . . , j, when the remaining budget available for the projects j + 1, j + 2, . . . , n is equal to t. Thus, we have:

F j (t) = max j h:=1 f h (x h ) s.t. j h:=1 x h ≤ A -t, x h ≥ 0, x h ∈ Z, h = 1, 2, . . . , j. (3)
We define F j (t) = 0 for t / ∈ [0, A], j = 1, 2, . . . , n and F 0 (t) = 0 for any t. Then we have the following recursive equations:

F j (t) = max x∈[0,A-t] {f j (x) + F j-1 (t + x)} = max 1≤k≤kj max x∈(t k-1 j ,t k j ] [0,A-t] {b k j -u k j t k-1 j + u k j • x + F j-1 (t + x)}, j = 1, 2, . . . , n. (4)
Lemma 1. All functions F j (t), j = 1, 2, . . . , n, are non-increasing on the interval [0, A].

The proof of this lemma immediately follows from definition (3) of the functions F j (t).

The DPA based on the equations ( 4) can be organized as follows. For each stage j = 1, . . . , n, for the state t = 0 we compute no more than

k j A values v x k = {b k j -u k j t k-1 j +u k j •x+F j-1 (t+x)}, 1 ≤ k ≤ k j , x ∈ (t k-1 j , t k j ] and put them into the corresponding lists L k . If we have for the next value v x k ≤ v x-1 k
, we do not put it into the corresponding list. This means that the values in each list L k are ordered in in non-decreasing order. For the next state t = 1, we only need to exclude the last element from the considered L k , k = 1, . . . , n, if it corresponds to an x which is not in the interval (t k-1 j , t k j ] [0, A -t] and compare the new k j last elements of the lists. If we continue in the same way for t = 2, . . . , A, we can calculate F j (t), t = 1, 2, . . . , A, in O(k j A) time. As a consequence the running time of the DPA using such a type of Bellman equations is O( k j A). A similar idea was presented in [START_REF] Shaw | An Algorithm for Single-Item Capacitated Economic Lot Sizing with Piecewise Linear Production Costs and General Holding Costs[END_REF].

The algorithms presented in [START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF] for the problem under consideration are based on the functional equations (3) and another technique to implement the graphical method. In contrast, the GrA presented in this paper is based on the equations (4).

Graphical algorithm

In this section, we develop an algorithm which constructs the functions F j (t), j = 1, 2, . . . , n, in a more effective way by using the idea of graphical approach proposed in [START_REF] Lazarev | A Graphical Realization of the Dynamic Programming Method for Solving NP-hard Problems[END_REF]. We will use the name DPA for the algorithm presented in Section 2 and GrA for this new algorithm.

The underlying general idea of improving the DPA for piecewise linear functions is as follows: Instead of going through all capacity values up to A, we keep the Bellmans functions depending on the capacity. They are again piecewise linear functions. Keeping these function pieces well organized for all capacities involves a considerable amount of technicalities and allows for running time improvements if the data representation is done in a clever way.

Below we prove that the functions F j (t), j = 1, 2, . . . , n, constructed in the GrA are piecewise linear. Any piecewise linear function ϕ(x) can be defined by three sets of numbers: a set of break points I (at each break point, a new linear segment of the piecewise linear function ends), a set of slopes U and a set of values of the function at the beginning of the intervals B. Let I[k] denote the k-th element of the ordered set I. The same notations will be used for the sets U and B as well. The notation ϕ.I[k] denotes the k-th element of the set I of the function ϕ(x). Then, for example, for x

∈ (t k-1 j , t k j ] = (f j .I[k -1], f j .I[k]], we have f j (x) = f j .B[k] + f j .U [k](x -f j .I[k]). Note that ϕ.I[k] < ϕ.I[k + 1], k = 1,
Ψ k j (t) = f j .B[k] -f j .U [k] • f j .I[k -1] + f j .U [k] • t + F j-1 (t)
according to Procedure 2.1.1.; 2.1.2. Construct the temporary function

Φ k j (t) = max x∈(fj .I[k-1],fj .I[k]] [0,A-t] {Ψ k j (t + x) -f j .U [k] • t} according to Procedure 2.1.2.; 2.1.3. IF k = 1 THEN F j (t) := Φ k j (t) ELSE F j (t) := max{F j (t), Φ k j (t)}. 2.
2. Modify the sets I, U, B of the function F j (t) according to Procedure 2.2.

3. The optimal objective function value is equal to F n (0).

The above algorithm uses Procedures 2.1.1. and 2.1.2. described below. In Procedure 2.1.1., we shift the function F j-1 (t) up by the value Given are k and j;

f j .B[k] - f j .U [k] • f j .I[k -1]
Ψ k j .I = ∅, Ψ k j .U = ∅ and Ψ k j .B = ∅.
FOR i := 1 TO |F j-1 .I| DO add the value F j-1 .I[i] to the set Ψ k j .I; add the value

f j .B[k] -f j .U [k] • f j .I[k -1] + f j .U [k] • F j-1 .I[i] + F j-1 .B[i] to the set Ψ k j .B; add the value f j .U [k] + F j-1 .U [i] to the set Ψ k j .U ;
Before describing Procedure 2.1.2., we present Procedure FindMax in which the maximum function ϕ(t) of two linear fragments ϕ 1 (t) and ϕ 2 (t) is constructed. The case when t * / ∈ Z is presented in Fig. 1. So, if both points t and t are integer, then ϕ.I contains only integer break points t. The running time of Procedure FindMax is constant.

Procedure FindMax

Given are the functions ϕ

1 (t) = b 1 + u 1 • t and ϕ 2 (t) = b 2 + u 2 • t
In the subsequent Procedure 2.1.2., we do the following. When we shift s to the right, we shift the interval I = [t lef t , t right ] of the length

f j .I[k]-f j .I[k -1].
We have to use the values Ψ k j (x) for x ∈ T to calculate Φ k j (t) at the point t = s . Since Ψ k j (x) is piecewise linear, it is only necessary to consider the values Ψ k j (x) at the break points belonging to T and at the end points of the interval T . So, if we shift s to the right by a small value x ∈ [0, ε] such that all the break points remains the same, then the value Φ k j (t) will be changed according to the value ϕ max (x).

:= Ψ k j .B[v] + Ψ k j .U [v] • (t lef t -Ψ k j .I[v -1]) -f j .U [k] • s ELSE b lef t := 0; 2.1.2.13. IF t right < A THEN b right := Ψ k j .B[w + 1] + Ψ k j .U [w + 1] • (t right -Ψ k j .I[w]) -f j .U [k] • s ELSE b right := 0;
{Ψ k j .B[s]+Ψ k j .U [s]•(Ψ k j .I[s]-Ψ k j .I[s-1])}-f j .U [k]•s ;
2.1.2.15. Denote function

ϕ lef t (x) := b lef t -(f j .U [k] -Ψ k j .U [v]) • x.
IF t lef t = A THEN ϕ lef t (x) := 0; 2.1.2.16. Denote function

ϕ right (x) := b right -(f j .U [k] -Ψ k j .U [w + 1]) • x. IF t right = A THEN ϕ right (x) := 0;
2.1.2.17. Denote function

ϕ inner (x) := b inner -f j .U [k] • x. IF T = ∅ THEN ϕ inner (x) := 0; 2.
F j .U [k], F j .U [k + 1] and F j .B[k], F j .B[k + 1], where F j .U [k] • (F j .U [k] -F j .U [k -1]) + F j .B[k] = F j .B[k + 1], (i.e., these
fragments are on the same line), then, to reduce the number of intervals |F j .I| and thus the running time of the algorithm, we can join these two intervals into one interval.

Procedure 2.2.

Given is F j (t);

FOR k := 1 TO |F j .I| -1 DO IF F j .U [k] = F j .U [k +1] AND F j .U [k]•(F j .U [k]-F j .U [k -1])+F j .B[k] = F j .B[k + 1] THEN F j .B[k + 1] := F j .B[k];
Delete the kth elements from F j .B, F j .U and F j .I; Lemma 2. Procedure 2.1.2. has a running time of O(|F j-1 .I|).

Proof.

Step [2.1.2.14] and the re-computation of T in step [2. 1.2.25] have to be performed with the use of a simple data structure. Let {q 1 , q 2 , . . . , q r } be a maximal subset of T having the following properties:

q 1 < q 2 < . . . < q r ;
there is no q ∈ T such that q i ≤ q < q i+1 and Ψ Analogously to Section 2, it is easy to show that F j (t), j = 1, 2, . . . , n, is a non-increasing function in t. Thus,

k j .B[q] + Ψ k j .U [q] • (Ψ k j .I[q] -Ψ k j .I[q -1]) ≥ Ψ k j .B[q i+1 ] + Ψ k j .U [q i+1 ] • (Ψ k j .I[q i+1 ] -Ψ k j .I[q i+1 -1]), i = 1, . . . , r -1. f 1 .I = {3, 10, 13, 25} f 2 .I = {5, 25} f 3 .I = {2, 4, 6, 25} f 4 .I = {3, 4, 25} f 1 .U = {0, 1, 1 3 , 0} f 2 .U = { 2 5 , 0} f 3 .U = {0, 2, 1 2 , 0} f 4 .U = {0, 0, 0} f 1 .B = {0, 0, 7, 8} f 2 .B = {0, 2} f 3 .B = {0, 0, 4, 5} f 4 .I = {0, 1, 4}
F j .B[k] ≥ F j .B[k + 1], j = 1, 2, . . . , n, k = 1, 2, . . . , |F j .I| -1.
Then, according to Procedure 2.2., there are no more than 2 • |F j .B| different values in the set F j .I.

Thus, the running time of the GrA is O k j min A, max j=1,2,...,n {|F j .B|} .

In fact, the running time is less than O( k j min{A, F * }), where F * is the optimal objective function value, since max j=1,2,...,n |F j .B| ≤ F * .

Example

Next, we will illustrate the idea of the GrA using the numerical example presented in Fig. 2. A full description of all calculations can be found in [START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF]. Here, we only present a short sketch. In this instance, we consider four projects with the profit functions f j (t), j = 1, 2, 3, 4 (see Table 1).

F 3 (t) = max{Φ 1 3 (t), Φ 2 3 (t), Φ 3 3 (t), Φ 4 3 (t)
F 3 .U = {0, -1 3 , -2 5 , -1 2 , -1, -2 5 , - 1 2
, -2 5 } and F 3 .B = {15, 14, 12, 11, 6 1 3 , 4, 1}. STEPS j = 4, k = 1, 2, 3 are performed in an easy way, i.e., to obtain the functions Φ 1 4 (t), Φ 2 4 (t) and Φ 3 4 (t), we have to shift the diagram of the function F 3 (t) to the left by the value 0, 3, 4 and up by the value 0, 1, 4, respectively. In Fig. 4, the maximum function F 4 (t) is presented.

To find an optimal solution at the point s = 0, we can do backtracking. We have x 4 = 4 and f 4 (x 4 ) = 4, x 3 = 6 and f 3 (x 3 ) = 5, x 2 = 5 and f 2 (x 2 ) = 2 as well as x 1 = 10 and f 1 (x 1 ) = 7. So, the optimal objective function value is F * (0) = 18.

In the GrA, we considered the following number of states s : 2+3+3+2 = 10 (for j = 1), 8 + 4 = 12 (for j = 2, where 4 states were considered for k = 2), 5 + 10 + 11 + 5 = 31 (for j = 3, where 5 states were considered for k = 1 and k = 4), 7 + 7 + 7 = 21 (for j = 4, i.e., during the shift of the diagram). So, in total we considered 10 + 12 + 31 + 21 = 74 states s . In the DPA, approximately 25(3 + 2 + 4 + 3) = 300 states will be considered. If we scale our instance to a large number M (i.e., we multiply all input data by M ), the running time of the DPA increases by the factor M , but the running time of the GrA remains the same. Of course, for each state in the GrA, we need more calculations than in the DPA. However, this number is constant and the GrA has a better running time.

An FPTAS based on the GrA

In this section, a fully polynomial-time approximation scheme (FPTAS) is derived based on the GrA presented in Section 3.

First, we recall some relevant definitions. For the optimization problem of minimizing a function F (π), a polynomial-time algorithm that finds a feasible 425 solution π such that F (π ) is at most ρ ≥ 1 times less than the optimal value F (π * ) is called a ρ-approximation algorithm; the value of ρ is called a worstcase ratio bound. If a problem admits a ρ-approximation algorithm, it is said to be approximable within a factor ρ. A family of ρ-approximation algorithms is called an FPTAS, if ρ = 1 + ε for any ε > 0 and the running time is polynomial polynomial. So, by using this technique, we can improve the FPTAS to have a running time of

O n • k j ε (1 + log log n) ,
Finally, we only note that an FPTAS based on a GrA was presented in [START_REF] Gafarov | A Graphical Approach for Solving Single Machine Scheduling Problems Approximately[END_REF] for some single machine scheduling problems.

Concluding Remarks

In this paper, we used a graphical approach to improve a known pseudopolynomial algorithm for a project investment problem and to derive an FPTAS with the best known running time.

The practical usefulness of the graphical approach is not limited to such a project investment problem or similar warehousing and lot sizing problems. The graphical approach can be applied to problems, for which a pseudo-polynomial algorithm exists and Boolean variables are used in the sense that yes/no decisions have to be made. This is the case for many applications of capital budgeting in science and engineering. However, e.g., for the knapsack problem, the graphical algorithm mostly reduces substantially the number of states to be considered but the time complexity of the algorithm remains pseudo-polynomial [START_REF] Lazarev | A Graphical Realization of the Dynamic Programming Method for Solving NP-hard Problems[END_REF]. On the other side, e.g., for the single machine scheduling problem of maximizing total tardiness, such a graphical algorithm improved the complexity from O(n p j ) to O(n 2 ) [START_REF] Gafarov | Transforming a Pseudo-Polynomial Algorithm for the Single Machine Total Tardiness problem into a Polynomial One[END_REF]. It seems to be challenging for future research to study other well-known combinatorial optimization problems with this approach to know exactly its performance for different problems and to deduce more general properties.

2 ,

 2 . . . , |ϕ.I| -1 and k j = |f j .I|. In each step j, j = 1, 2, . . . , n, of the subsequent algorithm, temporary piecewise linear functions Ψ i j and Φ i j are constructed. These functions are used to define functions F j (t), j = 1, 2, . . . , n, The functions F j (t) are piecewise linear as well.For t ∈ Z, their values are equal to the values of the functions F j (t) in the DPA. Let ϕ.I[0] = 0 and ϕ.I[|ϕ.I| + 1] = A. The points t ∈ ϕ.I and the other end points of the intervals with the piecewise linear functions considered in this article will be called break points. To construct a function in the GrA means to compute their sets I, U and B. Then the GrA is as follows. At each stage j = 1, . . . , n, we compute the temporary functions Ψ k j (t) and Φ k j (t) which are used to compute F j (t). The key features of the algorithm consist in the computation of Φ k j (t), modifying a piecewise linear function F j-1 (t) by changing linear fragments.Graphical algorithm1. Let F 0 (t) = 0, i.e., F 0 .I := {A}, F 0 .U := {0}, F 0 .B := {0}; 2. FOR j := 1 TO n DO 2.1. FOR k := 1 TO k j DO 2.1.1. Construct the temporary function

  and increase all slopes in its diagram by f j .U [k]. If all values t ∈ F j-1 .I are integer, then all values from the set Ψ i j .I are integer as well. It is obvious that Procedure 2.1.1. can be performed in O(|F j-1 .I|) time. Procedure 2.1.1.

Procedure

2. 1 . 2 .

 12 14. IF T = ∅ THEN b inner := 0 ELSE b inner := max s=v,v+1,...,w

STEP j = 1 ,

 1 k = 1. According to Procedure 2.1.1, we have Ψ k j (x) = 0, Ψ k j .I = {0}, Ψ k j .U = {0} and Ψ k j .B = {0}. Next, we consider each iteration of the cycle [2.1.2.5] in Procedure 2.1.2. First consider, s'=0. Before the first iteration, we have T = ∅, t lef t = 0,

Figure 2 :, Φ 3 1 and Φ 4 1 .

 21 Figure 2: Functions f j (t)

Figure 3 :Figure 4 :

 34 Figure 3: Calculations in the example

Figure 5 :

 5 Figure 5: Substitution of columns and modification of F l (t)

  2.1.2. IF w + 1 is not defined THEN let w + 1 = |Ψ k j .I|; 2.1.2.8. IF v is not defined THEN let v = |Ψ k j .I|; 2.1.2.9. IF t lef t < A THEN ε lef t := Ψ k j .I[v] -t lef t ELSE ε lef t := A -s ; 2.1.2.10. IF t right < A THEN ε right := Ψ k j .I[w +1]-t right ELSE ε right := +∞; 2.1.2.11. ε := min{ε lef t , ε right }; 2.1.2.12. IF t lef t < A THEN b lef t

	2.1.2.5. WHILE s ≤ A DO
	2.1.2.6. IF T = ∅ THEN let
	w + 1 = argmax i=1,2,...,|Ψ k j .I| {Ψ k j .I[i]|Ψ k j .I[i] > t right }
	and v = argmin i=1,2,...,|Ψ k j .I| {Ψ k j .I[i]|Ψ k j .I[i] > t lef t };
	2.1.2.7.
	2.1.2.1. Given are k, j and Ψ k j (t);
	2.1.2.2. Φ k j .I := ∅, Φ k j .U := ∅ and Φ k j .B := ∅;
	2.1.2.3. s := 0, t lef t := s + f j .I[k -1], t right := min{s + f j .I[k], A};
	2.1.2.4. Let T = {Ψ k j .I[v], Ψ k j .I[v + 1], . . . , Ψ k j .I[w]} be the maximal subset of
	Ψ k j .I, where t lef t < Ψ k j .I[v] < . . . < Ψ k j .I[w] < t right ,

Let T := {t lef t } T {t right };

  2.1.2.19. add the values from ϕ max .I increased by s to the set Φ k j .I; 2.1.2.20 add the values from ϕ max .B to the set Φ k j .B; 2.1.2.21. add the values from ϕ max .U to the set Φ k j .U ; 2.1.2.22. IF ε = ε lef t THEN exclude Ψ k j .I[v] from the set T and v := v+1; 2.1.2.23. IF ε = ε right THEN include Ψ k j .I[w + 1] to the set T and w := .1.2.25. t lef t := s + f j .I[k -1], t right := min{s + f j .I[k], A}, recompute T (for details, see the proof of Lemma 2); 2.1.2.26. Modify the function Φ k j according to Procedure 2.2. (described below). Next, we present Procedure 2.2. used in step [2.1.26] of Procedure 2.1.2. In Procedure 2.2., we combine two adjoining linear fragments that are in the same line. That means that, if we have two adjacent linear fragments which are described by the values (slopes)

	w + 1;
	2.1.2.24. s := s + ε.
	1.2.18. Construct the piecewise linear function
	ϕ max (x) := max x∈[0,ε] {ϕ lef t (x), ϕ right (x), ϕ inner (x)}
	according to Procedure FindMax;

2

Table 1 :

 1 Functions f j (t) Analogously to the proof of Lemma 4, after each step [2.1.3] of the GrA, the function F j (t), j = 1, 2, . . . , n, has only integer break points from the interval [0, A]. Each function Φ i j .I, j = 1, 2, . . . , n, i = 1, 2, . . . , k j , has only integer break points from [0, A] as well. So, to perform step [2.1.3], we need to perform Procedure FindMax on no more than A + 1 intervals. Thus, the running time of step [2.1.3] is O(A). According to Lemmas 1 and 2, the running time of steps [2.1.1] and [2.1.2] is O(F j .I), where F j .I ≤ A. The running time of step [2.2.] is O(F j .I) as well.

	Theorem 1. The GrA finds an optimal solution of the problem in
	O	k j min A,	max j=1,2,...,n	{|F j .B|}
	time.
	Proof.

Acknowledgments

This work has been supported by grants RFBR 13-01-12108, 13-08-13190, 15-07-03141, 15-07-07489 and DAAD A/1400328 and HSE Faculty of Economics.

We can keep track of the set {q 1 , q 2 , . . . , q r } by storing its elements in increasing order in a Queue Stack, i.e., a list with the property that elements at the beginning can only be deleted while at the end, elements can be deleted and added [START_REF] Aho | Data Structures and Algorithms[END_REF]. This data structure can easily be implemented such that each deletion and each addition requires a constant time. So, steps [2. 1.2.14] We remind that in the DPA, the functional equations (4) are considered. In fact, in Procedure 2.1.1., we construct the function

and in Procedure 2.1.2., we construct the function

Unlike the DPA, to construct Φ k j (t) in the GrA, we do not consider all integer

but only the break points from the interval, since only they influence the values of Φ k j (t) (and in addition t lef t , t right ).

Step [2.1.3] can be performed according to Procedure FindMax as well, i.e., to construct the function F j (t) := max{F j (t), Φ i j (t)}, their linear fragments have to be compared in each interval, organized by their break points. It is easy to see that we do the same operations with the integer points t as in the DPA. So, the values F j (t), t ∈ Z, are the same for the GrA and the DPA, and we can state the following: Lemma 3. The values F j (t), j = 1, 2, . . . , n, at the points t ∈ [0, A] Z are equal to the values of the functions F j (t) considered in the DPA. Lemma 4. All functions F j (t), j = 1, 2, . . . , n, are piecewise linear on the interval [0, A] with integer break points.

Proof. For F 0 (t), the lemma is true. In Procedure 2.1.1., all break points from the set Ψ i 1 .I are integer as well (see the comments after Procedure 2.1.1.). Since all points from f 1 .I are integer, we have ε ∈ Z and as a consequence, s ∈ Z.

According to the Procedure FindMax, all points ϕ max .I considered in Procedure 2.1.2. are integer. So, all break points from Φ i j .I, i = 1, 2, . . . , k j , are integer as well. Thus, the break points of the function F 1 (t) := max{F 1 (t), Φ i 1 (t)} are integer, if we use Procedure FindMax to compute the function max{F 1 (t), Φ i 1 (t)}. Analogously, we can prove that all break points of F 2 (t) are integer, etc.

It is obvious that all functions F j (t), j = 1, 2, . . . , n, constructed in the GrA are piecewise linear. Thus, the lemma is true. and b inner = -3 + 1 • (25 -0) -1 • 15 = 7. We get ϕ lef t (x) = 0 + (1 -1)x, ϕ right (x) = 0 and ϕ inner (x) = 7 -x. Thus, we obtain ϕ max (x) = 7 -x. We get s = s + 7 = 22; Finally, consider s'=22. We have T = ∅, t lef t = 25, t right = 22 + 10 = 32 and

and Ψ k j .B = {3 2 3 }. This step is performed analogously. We have to consider s = 0, 12, 15. We obtain Φ 3 1 .I = {12, 15, 25}, Φ 3 1 .U = {0, -1 3 , 0} and Φ 3 1 .B = {8, 7, 0}.

This step is performed analogously. We have to consider s = 0, 12.

We obtain Φ 4 1 .I = {12, 25}, Φ 4 1 .U = {0, 0} and Φ 4 1 .B = {8, 0}.

So, after STEP j = 1, we have

3 , -1, 0} and F 1 .B = {8, 8, 7, 0}, see Fig. 3.1. In fact, the function F 1 (t) is obtained only from the two functions Φ 2 1 and Φ 3 1 , where Φ 2 1 is the maximum function on the interval [0, 15] and Φ 3 1 is the maximum function on the interval [START_REF] Schemeleva | Lot-sizing on a Single Imperfect Machine: ILP Models and FPTAS Extensions[END_REF][START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF]. Next, we only present the states considered and the functions calculated in the steps j = 2, 3, 4. As mentioned before, the full description of all calculations can be found in [START_REF] Gafarov | A Graphical Approach to Solve an Investment Optimization Problem[END_REF]. Then we will be able to approximate the function F j (t) into a similar function with no more than 2 n 2 ε break points (see Fig. 5). Furthermore, for such a modified table representing a function F j (t), we will have

If we do the rounding and modification after each step [2.2], the cumulative error will be no more than nδ ≤ εF (π * ), and the total running time of the n runs of the step [2.2] will be

i.e., an FPTAS is obtained. In [START_REF] Chubanov | An FPTAS for a Single-Item Capacitated Economic Lot-Sizing Problem with Monotone Cost Structure[END_REF], a technique was proposed to improve the complexity of an approximation algorithm for optimization problems. This technique can be described as follows. Let there exist an FPTAS for a problem with a running time bounded by a polynomial P (L, 1 ε , U B LB ), where L is the input length of the problem instance and U B, LB are known upper and lower bounds, respectively. Let the value U B LB be not bounded by a constant. Then this technique enables us to find in P (L, log log U B LB ) time values U B 0 and LB 0 such that LB 0 ≤ F * ≤ U B 0 < 3LB 0 , i.e., U B0 LB0 is bounded by the constant 3. By using such values U B 0 and LB 0 , the running time of the FPTAS will be reduced to P (L, 1 ε ), where P is the same


Mohamed Kebdani
email: mohamed.kebdani@ec-lille.fr

Geneviève Dauphin-Tanguy
email: genevieve.dauphin-tanguy@ec-lille.fr.

Antoine Dazin
email: antoine.dazin@ensam.eu.

Patrick Dupont
email: patrick.dupont@ec-lille.fr.

Bond Graph Model of a mechanically Pumped Biphasic Loop (MPBL)

Keywords: Bond Graph, Biphasic, dynamic, Systems, Evaporator, cooling. I

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

INTRODUCTION

An electronic module implanted in embedded systems generates heat flows that require evacuation. For example the limit junction temperature of a silicon chip is 125 ° C; if it is exceeded a risk of malfunction is considered. Volume density of the power dissipated by Joule effect in a converter is quadrupled every four years, providing values that could reach 100 for the arrival of 2020. Serin, V. [START_REF] Serin | Hydrodynamic and thermal study of the vaporization in a square section of micro channel: application to two-phase micro capillary pumped loops[END_REF]. The cooling of electronic devices based only on exchanges by natural convection is designated to low power. Indeed, the apparent heat transfer coefficient remains low: [START_REF] Medjaher | Bond Graph Model of a Vertical U-Tube Steam Condenser Coupled with a Heat Exchanger[END_REF][START_REF] Mukherjee | Modeling and simulation of engineering systems through bond graphs[END_REF][START_REF] Muraoka | Experimental and theoretical investigation of a capillarity pumped loop with a porous element in the condenser[END_REF][START_REF] Bouamama | Supervision of an industrial steam generator. Part I: Bond graph modeling[END_REF][START_REF] Ould Bouamama | Bond Graph modelisation of steam condenser[END_REF][START_REF] Platell | Model Evaporator of two-phase fluid loop with capillary pumped[END_REF][START_REF] Schmidt | VDI tables constants steam[END_REF][START_REF] Serin | Hydrodynamic and thermal study of the vaporization in a square section of micro channel: application to two-phase micro capillary pumped loops[END_REF][START_REF] Thoma | Introduction to bond graphs and their applications[END_REF][START_REF] Wang | Capillary Heat Loop Technology: Space applications and recent Canadian activities[END_REF][13][14][15] . Furthermore, the finned radiators have apparent exchange coefficients up to 300 . For higher power, fluids loops are more adapted. These fluidic systems can be classified as follows:

i/ Mechanically Pumped Fluid Loop: MPFL. The loop is screwed on the circuit support. Inside the loop, a single phase refrigerant liquid is pushed by a pump. It is a forced convection whose apparent coefficient of heat exchange can reach . The fluid is used to extract heat at the electronic component and transport it to a cold source. ii/ Capillary Pumped Fluid CPF, which have been originally developed in the seventies, in the context of space applications. This system is based on the use of the latent heat of a fluid and is potentially able to evacuate important quantities of heat flows. The refrigerant movement is natural and obtained by capillary effect. CPF class is divided into two sub-categories developed in parallel. In the United States: CPL (Capillary Pumped Loop), developed by Stenger in 1966, Muraoka [START_REF] Muraoka | Experimental and theoretical investigation of a capillarity pumped loop with a porous element in the condenser[END_REF] in Russia, the LHP (Loop Heat Pipe), designed by Gerasimov and Maydanik 1975, Wang et al. [START_REF] Wang | Capillary Heat Loop Technology: Space applications and recent Canadian activities[END_REF]. The only difference between the two architectures is the positioning of the tank. CPF are distinguished by: transportation capacity powers ranging from 500W to 24000W, architectural flexibility and phases separation eliminating every phasic interactions.

iii/ Mechanically Pumped Biphasic Loop (MPBL).  MPBL components.

The principle of this system is quite similar to the previous one, but the fluid, trapped inside the loop, is moved by a pump. This allows greater flow rates and consequently greater heat flow transfer. The main components of MPBL are the following ones:

a. The evaporator. The heat source propagates heat flow, by conduction and convection, toward the fluid. This one extracts the flow by a phenomenon of vaporization, and evacuates it to the condenser. An evaporator is place of strongly and complex coupled thermal-hydraulics phenomena: evaporation, change of state, two-phase flow, dynamic movement of the evaporation front, evaporation rate...

b. The condenser.

Because of the heat exchange realized inside the condenser between the fluid source and a secondary refrigerant, the fluid condenses and releases the stored energy. Also, the condenser under cools the refrigerant to avoid cavitations.

c. The Tank and the pipes.

The thermo hydraulic control of the entire loop is possible through a thermal regulation of the fluid in the tank. This regulation allows the control of the pressure at the liquid / vapor interface in the tank. Knowing the losses, evaporation pressure is deduced. In general, the losses are small: it means that the pressure in the tank represents a reference for the whole system performance. In addition, the tank must be able to supply the loop with the liquid quantity needed for a good functioning, regardless whatever is the regime.

e. The pump.

The mechanical pump is located upstream of the evaporator (Fig. 6), the pump motor is powered by an external voltage source. Currently, we are in a step of feasibility study and performance evaluation of such a device.

II. BOND GRAPH MODEL. II.1. Generalized variables used in the model.

Two complementary pairs of variables, involved in the calculation of the transferred powers are defined in the sense of the bond graph methodology as:

For the hydraulic part, the effort variable is the pressure in the loop and the flow variable is the mass flow rate. The resulting model is a "pseudo" bond graph. The choice of the mass flow rate is justified by the fact that it is more adapted than volume flow rate for the formulation of the continuity equation in the case of a biphasic flow, in which the density is varying with time.

For the thermal part: the effort variable is the temperature and the flow variable is enthalpy flow. Enthalpy flow is preferred to entropy flow because energy balance is naturally written with this first variable. Furthermore, the choice of these variables allows constructing the thermodynamic functions of refrigerant used and calculating other system variables.

II.2. BG used to model a MPBL. II.2.1. Model assumptions.

The refrigerant used for the simulation is water; its thermo physical features are extracted from Schmidt [9]. The model is based on the following assumptions: The liquid flow is considered one-dimensional. Hydraulic inertia of fluid is not taken into account. Pressure is imposed in the upper part of the tank. The liquid-vapor mixture has uniform pressure. All pipes are perfectly adiabatic. There is no fluid storage in the pipes.

II.2.2. Word Bond Graph model of the loop.

The Word Bond Graph model is proposed Fig. 6 in which we find the loop components described in the introduction. 3) Sub-model 3: surface between the two previous regions. A multiport element R provides coupling between the two areas, and describes the phase change. Initially, the fluid is mono-phasic. After heat source activation, its temperature rises gradually to reach the saturation value. The fluid becomes biphasic. To specify this state change, and then be able to activate the model C of the evaporator, a Boolean indicator is introduced. This indicator noted ( ) is built according to the following logic: First, the saturation temperature corresponding to the pressure in the evaporator is calculated, by applying a polynomial approximation to the results provided by the thermodynamic tables: 11 If the liquid temperature is greater than , the fluid is biphasic and we take: otherwise .

 Hydraulic and heat balances inside the evaporator.

Conservation laws monitoring an evaporator operating in biphasic regime with phase change are:

Mass balance: 12 Energy balance : 13

Heat flow entering into the evaporator: 14

Heat flow leaving the evaporator: 15

Phase change equation: 16

 Determination of the mixture pressure " " The heat flow and the mass flow (Fig. 7) are calculated (eq. 5 and 7), and after integration (eq. 17 and 18) the specific enthalpy (eq.19) and specific volume (eq.20) are obtained. From these two values and the enthalpy diagram is raised the mixture pressure (eq.23). The mass of liquid in the evaporator is: 17

The enthalpy in the evaporator is: 18 From these values and the evaporator volume (

) we can obtain the specific enthalpy and volume: 19 20 Furthermore, the specific volume of the mixture is: 21 And the specific enthalpy of the mixture is: 22

After some arrangements we get the 12 th degree polynomial where the unknown variable is the pressure : All functions: are some known thermo-dynamics equations depending on pressure .

23

To calculate the pressure of the mixture, 20sim sends at each time, to Matlab the specific enthalpy h and the specific volume v. Using these two values Matlab solves the polynomial (eq. 23), and calculates the adequate root.

 Vapor quality "X".

The vapor quality expression is:

24

II.2.4. Model of the other elements of the MPBL.

At the present study, the pump is considered as a flow rate source, the tank is modeled as a source of pressure. Pipes /s. The temperatures in the evaporator (blue curve) and the condenser (orange curve) starts from the initial value 310K and increase very slightly, a few hundredths of a degree. This temperature increase is due to fluid friction.

III.2. Dynamic study of liquid phase thermal behavior.

In pure monophonic operating conditions, the MPBL is mainly a place of thermal process. In this section a dynamic study of the thermal behavior is proposed.

To perform an analytical study of temperatures in the exchangers, a simplified analogous electrical model of the MPBL (Fig. 10 (a)) and its BG (Fig. 10 (b)) are proposed, to enable better comprehension of the real MPBL behavior. Thermal flow sources are analogous to current sources Cold source" and (Fig. 10 (a)). The evaporator and condenser are reduced to capacitors noted successively: Cevap and Ccond. Pipes are assimilated to an electrical resistor. where:

cond liq p cond evap liq p evap m C C m C C _ _ . .   . 1 0 m p pipe q C R 
The different parameters are supposed to be constant and not temperature dependent. The model is linear. The differential state equation of the equivalent model is (eq.25):

                              u fr ch B cond evap x cond evap A cond cond evap evap x cond evap P P C C T T RC RC RC RC T T                                               1 0 0 1 1 1 1 1 With Initial condition:        cond evap T T x _ 0 _ 0 0 .
We can remark that:

, which means that the exchangers have an integrator-type behavior. The deduced transfer matrix is thus (eq.26): The stabilized regime is given by (eq.27 and 28):

                                                                 
              cond evap evap cond cond evap cond evap T T T T _ 0 _ 0 1 1 ) 1 1 ( 1 ) ( ) (     ) ( ) ( ) ( _ _ _ 0 _ _ 0 _

III.3. Cold and hot sources activated, no changing phase.

A second simulation, injected powers have opposite signs but same amplitude ( = ) is now considered. The temperatures are initially the same (293°K). A simple development of (Eq. 26) allows writing:

) 1 1 ( 1 ) 1 ( ) 1 1 ( ) ( ) ( _ 0 _ 0 0 cond evap cond evap evap cond cond evap evap ch evap s s T T s s C s P s T               ) 1 1 ( ) 1 ( 1 ) 1 1 ( ) ( ) ( _ 0 _ 0 0 cond evap cond evap evap cond cond evap cond ch cond s s T s T s C s P s T              
First of all, we can remark that integrator behavior has been cancelled. It is possible to predict the stabilized regime: In this case also, the graph 13 shows very similar temperatures evolution in both, electrical and MPBL model.

III.4. Phase change.

A final simulation with phases changes in the heat exchangers is now described and analyzed:

Analysis of the Figure 14:

The application of 2.0 kW (curve 1) at t = 1s leads to a rise of fluid temperature within the evaporator and the condenser (Curve 3, 4). The evaporator temperature reaches the saturation value of water T = 373K at t = 7s. Curves 8 and 9 show the evolution of the mixture temperature in both evaporator and condenser. The curve 5 shows the state indicator that displays the value 1 at the instant when the liquid temperature is: 373K. By switching from "0" to "1", it informs about the moment of saturation, namely t = 7s. The fluid is now biphasic: liquid/vapor mixture, whose vapor quality is up to 0.3 %. Curve 10 (in evaporator) and curve 11 (in condenser). The curve 6 and 7 trace the evolution of the mixture pressure respectively in the evaporator and condenser. It is noted that once the cooling triggered at t = 15s (curve 2) the steam quality (curves10 and 11) and the pressure (curve 6 and 7), start decreasing before stabilizing. The steam quality converges to an equilibrium value 0.2 %. Note that with a power = 2000 W and a flow rate = 5 /s, the specific enthalpy can be calculated as given: kJ / Kg. According to Mollier diagram (P, H); theoretical corresponding steam quality would be 2%, a value which coincides well with the results available at the model output.



  Existing Numerical models. Few papers in the literature deal with MPBL components bond graph (BG) modeling. A BG model detailed of a condenser is proposed by Medjaher. and al. [3]. The work is about the functioning of the heat exchanger, and presents a physical model formulation of the problem with a BG whose simulation results are consistent with the experimental. The work of B. Bouamama and al. [7] fits along the same lines of modeling a steam condenser with BG. The scientific literature shows few studies that provide numerical models of fluid loops, and even less of two-phase loops. However, we can cite the work of B. Ould Bouamama and al. [6] studying the complex behavior of a steam generator consisting of a condenser, steam accumulator and pipes. The authors bring forward a comprehensive model BG of the installation. Another numerical modeling work is proposed by Platell and al. [8] where they present a non-stationary model of CPL. Also in the same vein, Lachasagne [2] proposes an analytical model dedicated to CPL and based on an original modeling approach called nodal. As to MPBL, at our knowledge, no modeling work of the complete loop has been published at present time. The work presented in this paper is intended to contribute to the theoretical study of a MPBL with a specific attention on the phase changes in the heat exchanger. Modeling and simulation are performed using the software 20sim by adopting the BG methodology (Thoma [11], and Mukherjee [4]). It is used to show up the power exchanges in a multidisciplinary dynamic system (mechanical, electrical, hydraulic, thermal ...) based on the principle of power conservation and causality property [1]. The control of mass and enthalpy balance is capital and must be completed for each component of the loop. The notion of junction "0 and its dual 1" facilitates the representation of mass and enthalpy conservation. The purpose is to be able to predict the behavior of the loop in transient and steady operations and consequently to provide a pertinent tool for the design of such loops. The results of the model are also used to analyze the physical behavior of the system during several operating conditions.

Fig. 6 :

 6 Fig.6: Word Bond Graph of a MPBL. II.2.3. Detailed Bond Graph of the evaporator. The global model Fig 7 of the evaporator is composed of three sub models: 1) Sub-model 1: Purely Liquid region of the evaporator. The model used is a passive element multiport RC, where both hydro and thermal energies are coupled: R describes the energy dissipation due to load losses. C expresses the heat storage in the thermal part. 2) Sub-model 2: two-phase mixing region. The model used is a passive two-port element C, storing mass in the hydraulic part and heat in the thermal part.

Fig. 7 : 3 Case 2 : 8 Initial enthalpy : 9 part : 10 

 7328910 Fig.7: BG Model of the evaporator.a. Energy equationThe general expression of the enthalpy flow is:

1 .

 1 modelise load losses (elements R) and fluid heating caused by frictions forces (elements S). At last, the condenser model is quite similar to the evaporator model (Fig 8). III. RESULTS AND ANALYSIS. Geometrical dimensions of MPBL components. III.1. Cold and hot thermal power sources disabled. Simulation conditions: Injected power = 0 W; Dissipated power = 0W. Initial pressure =1.3 bar; initial temperature = 310K Mass flow rate = 5

Fig 9 :

 9 Fig 9: Temperature evolution in 1) Evaporator (orange), 2) Condenser (bleu).

Fig. 10 .

 10 Fig.10. (a) Equivalent electrical model of the MPBL, (b) BG of the equivalent MPBL model.



  Analysis case n°1: Hot source = 0W, Cold source = 0W. The initial temperatures are different (293°K /313°K).

  The permanent regime behavior depends only on liquid masses and initial temperatures in exchangers. Analysis case n°2: Hot source = 200W, Cold source = 0W. The temperatures are initially the same (293 °K). According to the system equation (Eq. 26) we have:

Fig. 11 .

 11 Fig.11. Transfer functions of the equivalent model.

Fig. 12 .

 12 Fig.12. Temperature evolution of both electrical model and MPBL. The graph 12 shows very similar temperatures evolution in both, electrical and MPBL model. The difference in the transient evolution of and is due to the influence of the filter Gcond(s) whose expression is shown in Fig.11.

  This means that final temperature in the evaporator is different from the one in the condenser (fig.13).

Fig. 13 .

 13 Fig.13. Applied power. MPBL Model. Equivalent Electrical Model

Fig. 14 :

 14 Fig. 14: Phase change in evaporator (red) and Condenser (Blue).

IV. CONCLUSION

The functioning principle of the cooling loops is based on the use of physical and thermodynamic properties of fluids. The concept of phase change and its efficiency in transferring energy was discussed. A second part is devoted to the presentation of a biphasic fluid loop model, in which the fluid is mechanically pumped. The modeling is based on Bond Graph approach. The third part includes the cases of simulation launched with the developed model. Note that this model is configurable, which means that it is possible to change, at any moment, the geometric conditions of the loop as well as the fluid properties. A test bench of the entire loop is under development, the purpose is validating the model in the short term with experimental results.

Variables