Robin Smith
email: robin.smith@manchester.ac.uk

Megan Jobson

Lu Chen

Recent Development in the Retrofit of Heat Exchanger Networks

Keywords: Heat exchanger network, Network pinch approach, Mixed-integer non-linear optimisation, Retrofit design

This work presents a methodology for heat exchanger network (HEN) retrofit, which is applicable to complex industrial revamps, considering existing networks and constraining the number of modifications.

The network pinch approach [1] has been modified and extended to apply to the HEN design in which the thermal properties of streams are temperature-dependent. The modified network pinch approach combines structural modifications and cost optimisation in a single step to avoid missing cost-effective design solutions.

Introduction and previous work

Numerous investigations have been carried out to improve the performance of HENs. The HEN design approaches can be grouped into three major categories: pinch analysis methods, mathematical programming methods and stochastic optimisation methods. Linnhoff and Hindmarsh [START_REF] Linnhoff | The pinch design method for heat exchanger network[END_REF] developed a method to design for the minimum hot and cold utility demand, for a set of hot and cold streams and a selected minimum temperature approach (ΔT min ), based on several heuristic rules. Because pinch analysis is difficult to apply to large scale problems, mathematical modelling methods were developed.

Optimisation of HEN is generally formulated as a mixed-integer non-linear programming (MINLP). To avoid being trapped into local optimum, simplifications are made to convert the MINLP problem into Linear programming (LP) problem [START_REF] Papoulias | A structural optimisation approach in process synthesis -I[END_REF], or Non-linear programming (NLP) problem as in the work of Yee et al. [START_REF] Yee | Simultaneous optimisation models for heat integration -I. Area and energy targeting and modeling of multi-stream exchangers[END_REF]. Sorsak and Kravanja [START_REF] Sorsak | MINLP retrofit of heat exchanger networks comprising different exchanger types[END_REF] and Ma et al. [START_REF] Ma | Constant approach temperature model for HEN retrofit[END_REF] have formulated MINLP models for the retrofit design of HEN. Compared to deterministic methods, using stochastic methods will give more chances to find the global optimum for MINLP problems, due to the random nature of the optimisation methods.

Commonly used algorithms in the design of HENs are Genetic Algorithms [START_REF] Ravagnani | Heat exchanger network synthesis and optimisation using genetic algorithm[END_REF] and Simulated Annealing [START_REF] Athier | Synthesis of heat-exchanger network by simulated annealing and NLP procedures[END_REF].

Non-constant thermal properties often arise when multi-component streams are cooled down or heated up, such as in refining preheat trains. Only few methodologies considered the varying thermal property (e.g. heat capacity) issue of process streams [START_REF] Grossmann | A rigorous disjunctive optimization model for simultaneous flowsheet optimization and heat integration[END_REF][START_REF] Castier | Energy targeting in heat exchanger network synthesis using rigorous physical property calculations[END_REF][START_REF] Ponce-Ortega | Optimal synthesis of heat exchanger networks involving isothermal process streams[END_REF]. These proposed methodologies have several limitations. In the work of Grossmann et al. [START_REF] Grossmann | A rigorous disjunctive optimization model for simultaneous flowsheet optimization and heat integration[END_REF], the network configuration has been fixed and not optimised. Castier and Queiroz [START_REF] Castier | Energy targeting in heat exchanger network synthesis using rigorous physical property calculations[END_REF] only estimated the minimum energy requirement, and detailed HEN designs were not

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 4
considered. The MINLP model proposed in the work of Ponce-Ortega et al. [START_REF] Ponce-Ortega | Optimal synthesis of heat exchanger networks involving isothermal process streams[END_REF] is for HEN synthesis only and not applicable to retrofit.

The network pinch approach [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF] combines physical insights into retrofit problems and mathematical programming techniques. The bottleneck of the existing network configuration is first identified by redistributing the heat loads of existing exchangers, which is referred to as pinching the existing network.

Then each candidate structural modification that may overcome the bottleneck of the HEN configuration is optimised at a time for maximum heat recovery. A list is generated after all suggested modifications are optimised, showing the corresponding maximum heat recovery for a given modified HEN topology. The difficult MINLP problem is then decomposed into Mixed-integer linear programming (MILP) problem and NLP problem. Although it is a sequential approach, it explores possible topology modifications in a systematic way and at the same time allows user interactions in the design procedure. This characteristic makes the network pinch approach to be a promising retrofit design methodology in industrial practice.

However, the existing network pinch approach assumed constant thermal properties (e.g. heat capacity)

and stream split fractions are not considered in pinching existing networks. Moreover, the existing approach only carries out cost-optimisation after the diagnosis stage. The design with minimum cost cannot be guaranteed since the selection of the potential modifications is not based on costs but energy demands. In this paper, the network pinch design method [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF] is modified to overcome these limitations.

Structure representation and modelling of heat exchanger network

There are some streams for which the thermal properties (e.g. heat capacity) are highly dependent on temperature. For those streams, multi-segment formulations are employed. The whole temperature range

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 5
is broken into several intervals (the stream in each interval is referred to as a segment) in which the thermal properties are assumed constant.

The representation of HEN structure in this work employs the node-based representation. As shown in Figure 1, the links between every component in the heat exchanger network are represented by nodes.

There are four nodes assigned to a heat exchanger unit (for example, process exchanger and utility exchanger): hot side inlet node, hot side outlet node, cold side inlet node and cold side outlet node. For unit operations, as they are only relative to one stream, two nodes associate with them, namely inlet and outlet nodes. Unit operations are devices designed to alter either the temperature or heat content of a stream. They represent the background process of the HEN without considering the actual details of the behaviour of the unit operation. For example, in preheat trains, the temperature of crude oil decreases after passing through a desalter. The desalter can be modelled as a unit operation. For a unit operation, temperature change of the stream is specified. Each node is associated with a unique temperature, which means that a new node is defined only if the temperature varies (e.g. the temperature of a stream branch after a splitting is the same as that before splitting). Therefore there is one node associated with supply temperature of each stream and also one node is associated with each stream splitter. The first node of each stream is always associated with the stream supply temperature; and the last node of each stream is always associated with the stream target temperature. Owing to the unique one-node-one-temperature data structure, redundant data are avoided and very complex networks are represented using moderate memory: stream splitting is considered in the HEN simulation and design; utility exchangers can be located at any places in the HEN; the number of exchangers on a given stream or branch thereof is not be inherently limited.

The HEN is modelled as an interconnected set of network elements, namely process heat exchangers, utility heat exchangers, stream splitters and mixers and unit operations. Because multi-segmented stream

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6
data are employed, the equations are not linear with respect to temperature. Therefore, with specified heat load of each exchanger unit, the node temperatures in the network are calculated sequentially. The calculation order of the node temperature is from the start of each stream to the end of that stream. After all the temperature nodes are calculated, the required heat transfer area is calculated for each exchanger.

The needed area is then compared to the existing area of each exchanger. Capital cost will be incurred if the needed area is larger than the existing one.

Network pinch approach for HEN retrofit

In the first step of the modified network pinch method, the network pinch is identified by redistributing heat transfer loads between the existing matches and in the meantime varying stream split fractions for maximum heat recovery (minimum utility demand). Both heat loads and stream split fractions are adjusted to make sure that the network pinch is not caused by the heat transfer area limits but the topology of existing HEN. The redistribution of heat loads and variation of stream split fractions in the network for maximum heat recovery is of great importance due to its effect on determining the bottleneck of existing HEN topology and suggesting structural modifications that may overcome the bottleneck. Given a network having N HX heat exchangers and N BR stream branches, the objective in the pinching problem is to minimise utility demand of the network (Utility), by varying heat loads of existing exchangers

( HX N Q Q Q  , , 2 1
) and stream split fractions (

BR N ff ff ff  , , 2 1
), subjected to minimum temperature approach constraints and stream enthalpy constraints (Equation 1).

  0 0 0 : , , , , , min , min min 2 1 2 1            j j calc out l in l in l out l N N TT TT T TH TC T TH TC st ff ff ff Q Q Q f Utility BR HX   (1) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

The number of variables N

VARY in the pinching problem is the sum of N HX (the number of heat exchangers) and N BR (the number of branches). Considering the large number of optimisation variables and the nonlinearity of the problem, the pinching problem is very complex. Normally N HX is significantly larger than N BR , since N BR is limited for the reason of controllability. The relative importance of the heat exchangers suggests decomposing the pinching problem into two levels (Figure 2). In the outer loop the stream split fractions are optimised, and in the inner loop the heat loads of existing exchangers are redistributed. The minimum utility requirement obtained in the inner loop serves as the objective function value for the outer loop optimisation. By implementing the two-level pinching procedure, the size and complexity of the pinching problem are reduced. As will be explained in following sections, the inner loop optimisation is formulated as a least-squares problem, which has an easier convergence and less calculation iterations compared with other types of NLP problems [START_REF] Nocedal | Numerical Optimisization[END_REF]. Moreover, solving least-squares problems will less depend on the quality of the initial guess, which is not easy to be always sensible provided in a complex NLP problem. As the inner loop of the pinching problem contributes more than the outer loop, the penalties incurred by decomposing the pinching problem on the optimality are small. In a decomposed way, there are more chances to find a better solution, in terms of robustness of the optimisation method due to less dependence on the initial guess and easier convergence.

The successive quadratic programming (SQP) method is a popular approach for process optimisation, and are most efficient if the number of active constraints is nearly as large as the number of variables, that is, if the number of free variables is relatively small [START_REF] Nocedal | Numerical Optimisization[END_REF]. The SQP algorithm is selected as the optimisation solver for the outer loop where only the flow rates of stream branches are optimised. In particular, the In the inner loop, the flow ratios of each branch respect to the main streams are fixed, thus the optimisation problem is related purely to heat loads, and subject to minimum temperature approach constraints and stream enthalpy balance constraints. The inner loop optimisation is a highly constrained non-linear programming problem, involving a large number of optimisation variables. The commonly used SQP algorithm, which is normally suit optimisation problems with small number of variables, is not applied in the inner loop, where a large number of variables, HX N , are involved. This work presents a novel approach for solving the inner loop optimisation problem.

As in this work, the limitation of assuming constant heat capacity with temperature is overcome. In the inner loop, the heat loads of existing exchangers are optimised toward minimum energy requirement. In order to solve this optimisation problem, first of all, a correlation needs to develop to associate the node temperature ( i T ) with heat loads of exchangers ( Q ) and the supply temperature ( j TS ) of the stream where the node i is located. A polynomial correlation is proposed to represent each node temperature in terms of enthalpy change (Equation 2).

j i j i j i j i j i E DH D DH C DH B DH A T          2 3 4
(2) where i DH is the accumulated enthalpy change from the start of the stream ( j TS ) to the particular node i, as shown in Equation 3:

   HX N l lj l il i ff Q y DH 1 (3) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 9 where 1 or 0  il y
, denotes the existence of exchanger unit l on the section of the stream from node i to the start of stream j; l Q represents the heat load of exchanger unit l; lj ff represents the flow ratio of the branch to the main stream j where the exchanger unit l is located.

Although there are other types of correlations that are commonly used (e.g. exponential, moving average), the polynomial correlation is proposed based on the following considerations:

 The formulation of polynomial correlation is flexible. The order of the formation can be adjusted to best suit the set of data. In this work, fourth order is proposed based on experiments, which gives the best fitting.

 The formulation of the polynomial correlation is simple. The analytical form of the first and second order derivatives can be attained easily. This allows employment of deterministic methods, which require detailed formulation and derivatives of the problem, in solving the inner loop optimisation problem.

 In the simulation of HEN, no formulation between the temperature and thermal properties are required since the stream data with varying thermal properties are input as multiple linear segments.

It is simple to attain the coefficients in Equation 2 from the multi-segmented stream data due to the popularity of the polynomial correlation and availability of the regression in many commercial tools.

The proposed polynomial correlation is implemented to optimise the heat loads in the inner loop. If there are utility heat exchanger units located on stream j, the overall cumulative enthalpy change applied in Equation 2 excludes heat loads of those utility heat exchanger units. Then the calculated temperature of stream before utility exchangers (TT ppcalc,j ) is equivalent to the temperature after process-to-process heat recovery. It is clear that the smaller the difference between TT ppcalc,j and the target temperature TT j , the lower the utility demand. Therefore, the utility consumption in terms of heat loads of process exchangers

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

10 can be formulated by Equation 4. In the equation, the minimum temperature approach constraints and stream enthalpy balance constraints are considered by associating penalty functions in the objective.

      } 0 , , min min{ 1 2 , 1 2 min min 1 2 ,                            SUT ST

TT

, is calculated by Equations 2, so each of the term in Equation 4 is a function of heat loads of all the process exchangers;  is the penalty factor for violated constraints.

It can be seen that by employing the polynomial correlation, the optimisation in the inner loop of the decomposed network pinching is formulated as a least-squares problem (Equation 4), which has an easier convergence and less calculation iterations compared with other types of NLP problems [START_REF] Nocedal | Numerical Optimisization[END_REF]. Levenberg-Marquardt algorithm [16] is used to solve the inner loop optimisation where heat loads of existing matches are distributed for maximum heat recovery.

The Levenberg-Marquardt algorithm is explained here in short. The searching process starts with an initial guess of 0

x , then  

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 11

In the work of Asante and Zhu [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF], each modification is modelled and optimised for maximum energy recovery and ranked according to the utility consumption. Up to this step, the diagnosis stage is accomplished, and followed by the optimisation stage which optimises the selected candidate modifications for minimum total annualised cost. However, in the diagnosis stage, capital cost is not considered, the heat loads may be redistributed such that additional area is needed in many existing heat exchanger units. In practice, adding additional heat transfer area to one heat exchanger unit is not only about installing new area. The associated pipe work also requires modifications, which is normally more expensive than the installed area [18]. For example, adding 500 m 2 to 5 existing units is much more expensive and less favourable in industry than adding same amount of 500 m 2 to a single unit.

To take into consideration retrofit costs, in the current study, the network pinch approach of Asante and Zhu [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF] is modified as shown in Figure 3. In the new network pinch approach, the structure searches in the diagnosis stage and the cost-optimisation are combined into one stage. The design problem then becomes a search for the most cost-effective structural changes in only one step, rather than sequential steps of a search for structural changes followed by a capital-energy optimisation, as employed in the original network pinch approach of Asante and Zhu [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF]. The one-step design approach avoids missing potentially cost-effective designs in the diagnosis stage, by ranking the alternative designs based on costs, rather than heat recoveries.

In the proposed one-step approach, each candidate modification is optimised for minimum total cost directly. The minimum total cost found in the optimisation of each modification is used later as the criterion for ranking potential modifications. The design problem thus becomes a non-linear programming  Compared with using deterministic methods, using stochastic methods will give more chances of finding global optimum for non-linear problems due to their random characteristics. The nonlinearity of the problem is increased greatly by the implementation of more accurate non-constant thermal properties for process streams.

 The network feasibility may be violated when simulating candidate solutions. It is not suitable to use deterministic methods, where a feasible initial design is normally required. The SA algorithm proposed in the work of Chen et al. [19] is capable of dealing with an infeasible initial design facilitated by the feasibility solver.

There are many advantages to employ the SA algorithm as the optimisation method, such as solutions of better quality, more constraints that can be handled, and fewer simplifications to the models of HENs to ensure convergence. However, relative longer computation times are the cost of these benefits. Compared to the original sequential order of structure searches in the diagnosis stage and the cost optimisation stage where a deterministic method is employed to solve the optimisation problem, more time is needed (at least 102 more time is required, the length of calculation time depends on the setting of SA parameters, see Appendix B.2) using the SA algorithm in minimising the total annualised cost of modifications that may overcome the bottleneck of the exiting network.

In order to reduce the calculation time, the SA algorithm [19] is modified such that only one run is needed in searching one type of structural modification. In retrofit design of HEN, there are four types of

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

13 structural modifications to overcome the network pinch by moving heat from below the pinch to above the pinch of the existing network. The four possible types of changes are: re-piping existing matches, resequencing existing matches, inserting a new match and introducing additional stream splitting to the existing network. In the optimisation, both continuous variables and structural options of a particular type are optimised. If the SA algorithm makes a structural move, the structural modification is selected from those suggested after pinching the existing network. In the SA optimisation procedure, the best design of each potential modification is stored and reported at the end of the structure searching process. Note that the best design for each structural change from the modified SA algorithm is slightly worse than those gained from running the SA algorithm hundreds of times to optimise each potential modification and identify the most cost-effective design with the particular modification. However, given N candidates in a type of structural modifications, the calculation time is reduced to only 1/N of that required in running the SA algorithm each time in optimising each single candidate.

Normally, the objective in HEN design is to design a practical cost-effective network. Therefore, In SA algorithm, total annualised costs, comprising utility costs and annualised capital costs, is taken as the objective function in current work.

The modified network pinch approach takes into account different levels of complexity in the HEN retrofit. The complexities addressed include non-constant thermal properties, the effects of stream split fractions on determining the network pinch and the combination of the diagnosis stage and the cost optimisation stage to avoid missing cost-effective designs. In order to reduce the computational time, different methods are used to solve different types of design problems. For instance, if the constant heat capacity assumption stands, and initially there are no stream splits in the existing network, the pinching problem is solved by a linear solver (e.g. E04MBE of the NAG FORTRAN library) rather than the SQP method which is particularly suitable for non-linear system; Figure 4 shows the overall structure of the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

14 modified network pinch approach together with Figure 5 showing the sub-section of the overall network pinch approach.

Case study

The application of the modified network pinch is illustrated through adding a new exchanger unit in an existing preheat train. The crude oil preheat train is the heat exchanger network that is associated with an existing atmospheric crude unit [20]. The multi-segmented stream data considering varying heat capacities (target temperature, supply temperature, enthalpy change) are shown in Table 1. The existing HEN structure is shown in Figure 6, while the existing heat transfer areas for exchanger units are presented in Table 2. The utility costs are calculated using unit cost data (annual cost per unit of energy, $/kW•y) and the calculated demand for each utility. The unit cost of hot utility is calculated from fuel cost, fuel net calorific value and efficiency of the atmospheric furnace. The unit cost of cold utility is taken from Gadalla [20]. The capital cost of heat exchanger modifications are calculated from simple correlations, which link a cost per unit area [18]. Unit cost of utilities and correlations to calculate heat exchanger capital cost from areas are shown in Table 3. The structural bottleneck of the existing network for ΔT min of 30 °C is identified using the developed two-level pinching method. The identified heat exchanger units at pinch location of the original network are shown in Figure 7. The change in flow fractions of braches, relative to the main stream, and redistribution of heat loads of existing matches for the pinched network are presented in Tables 5 and6, respectively. The existing crude oil preheat train is also pinched using the SQP algorithm to optimise the heat loads and stream split fractions in a single step.

In addition, the network is pinched by using the method of Asante and Zhu [START_REF] Asante | An automated approach for heat exchanger network retrofit featuring minimal topology modifications[END_REF], in which only heat loads are redistributed for maximum heat recovery.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 15

The minimum utility demands attained using the three methods are compared in Table 4. Table 7 shows the performance of the optimised design with a new match using the developed approach. The performance of the best design without any topology modifications is also shown in Table 7, and compared with the results of adding a new exchanger.

Table 4 shows that the developed two-level approach gives pinched HEN with most heat recovery, indicating that it is more capable of exploring the scope for recovering heat with the existing HEN configuration and determining the network topology bottleneck. It can be seen from Table 7 that by adding a new exchanger, the cost-effective retrofit design saves 20,358 kW hot and cold utilities, around 23% of the base case energy consumption. Table 7 also shows that a further 4,376 kW hot utility can be saved by adding a new exchanger, compared with no structural changes. These results indicate that in some cases, cost-effective retrofit does not always require topology modifications. Figure 8 presents the most cost-effective design with a new exchanger.

Conclusions

A retrofit design methodology is proposed for HENs of process streams with temperature-dependent thermal properties. The two-level pinching approach is developed for the optimisation of all continuous variables so that the heat recovery of the existing HEN is exploited to make sure that the bottleneck is the network topology rather than heat transfer areas. Moreover, the search for structural changes and capitalenergy optimisation are combined into a single step for the first time, in order not to miss cost-effective designs. The new HEN retrofit design approach enables identifying the most critical HEN configuration changes, and provides access to designers, which makes sure the retrofit designs are mathematically optimum and industrially applicable. [13] Numeric Algorithm Group, NAG FORTRAN library, Routine E04UCF, Oxford, UK, 1990.

[14] J.E. Dennis, R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice-Hall Inc., New Jersey, USA, 1983.

[15] L.T. Biegler, I.E. Grossmann, A.W. Westerberg, Systematic methods of Chemical Process Design, Prentice Hall Inc., USA, 1997

[16] P.E. Gill, W. Murray, J. SIAM, Numer. Anal., 15(5) (1978) 977-992.

[17] D. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, SIAM Journal on Applied Mathematics, 11 (1963) 431-331.

[18] R. Smith, Chemical process design and integration, John Wiley & Sons Ltd., 2nd Edition, West Sussex, England, 2005.

[19] L. Chen, M. Jobson, R. Smith, Automated heat exchanger network design for process streams with temperature-dependent thermal properties, AIChE spring meeting, paper 165f, New Orleans, 2008.

[20] M.A. Gadalla, Retrofit of heat-integrated crude oil distillation systems, PhD Thesis, UMIST, Manchester, UK, 2002

AQT 10 

 10 Total existing area of an existing network (m 2 ), in Figure3cu Cold utility i DH Accumulated enthalpy change from the start of the stream to the particular node i (kW), in Equations 2, defined by Equation3exist E Energy demand of an existing network (kW), in Figures 3 lj ff Flow ratio of the branch to the main stream j where the exchanger unit l is located on HEN Heat exchanger network hu Hot utility LP Linear programming MILP Mixed integer linear programming MINLP Mixed integer non-linear programming mods Modifications, in Figures 3 BR N Number of stream branches in a heat exchanger network HX N Number of heat exchangers in a heat exchanger network ST N Number of process streams in a heat exchanger network SUT N Number of process streams on which there is a utility heat exchanger, in a heat exchanger network y N var Number of variables in the pinching problem NLP Non-linear programming l Heat load of heat exchanger l (kW) i R m ax Maximum heat recovery for the network with i modifications (kW) in Figures 3 SA Simulated annealing SQP Successive quadratic programming i Temperature of node i (°C)in l TC Inlet temperature of cold stream in heat exchanger l (°C) out l TC Outlet temperature of cold stream in heat exchanger l (°C) in l TH Inlet temperature of hot stream in heat exchanger l (°C) out l TH Outlet temperature of hot stream in heat exchanger l (°C) j TS Supply temperature of stream j (°C) j TT Target temperature of stream j (°C) j calc TT , Calculated target temperature of stream j, including transfer heat with utility (°C) j ppcalc TT , Calculated temperature of stream j after process to process heat transfer in a heat exchanger network (°C) il y Binary variable, denotes the existence of exchanger unit l on the section of the stream from node i to the start of stream Damping factor in Levenberg-Marquardt algorithm

  are the hot inlet temperature and cold inlet temperatures of heat exchanger l ; outlet temperature and cold outlet temperatures of heat exchanger l ; j calc TT , is the calculated target temperature of stream j, and j TT is the specified target temperature of stream j.

  M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8 optimisation subroutine E04UCF of the NAG FORTRAN library [13] is used to implement the Quasi-Newton method [14] in the solution. The details of the SQP algorithm are explained in Biegler et al. [15].

  is calculated by Equation 4. The Jacobian matrix   0 x J is then calculated. The searching direction and length of the step x  is determined by solving   until the final solution is found. The setting of initial damping factor 0  and  are more or less experimental and problem specific [17]. In this work, these values are taken as 0.001 and 10.0 respectively based on trial and error. After the bottleneck of the HEN configuration is diagnosed by pinching the network, different types of modifications are suggested to overcome the pinch. The modifications fall into four groups: re-piping existing matches, re-sequencing existing matches, adding new matches and introducing stream splitting.

( 12 

 12 NLP) problem. The simulated annealing algorithm with a feasibility solver [19] is employed to solve this NLP problem because of the following considerations: In the SA algorithm, the simulation and the optimisation algorithm are decoupled. No derivatives are required. The NLP optimisation problem is treated as a black-box; only objective function values of trial solutions are needed. That is, no simplifications to the simulation models are required for the sake of optimisation. The calculation of capital cost and operating cost can also be as complex as is required.

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 8

 1234568 Figure 1 HEN structure representation

Figure8

  Figure1 HEN structure representation

Table 1

 1 Segmented stream data of illustrative example 4.1

			Supply	Target
	Stream	Name	temperature	temperature	Enthalpy change
			°C	°C	MW
	1	Pump-Ar 1	298	268	12.828
	2	Bott Cool 1	339	299	9.604
			299	259	9.164
			259	219	8.705
			219	179	8.228
			179	139	7.735
			139	100	7.024
	3	Pump-Ar 2	250	210	14.416
			210	200	3.467
	4	Bott Cool 2	257	217	1.142
			217	177	1.077
			177	137	1.011
			137	97	0.943
			97	57	0.874
			57	50	0.148

Approach 1: Network pinch: only heat loads varied; Approach 2: SQP algorithm: heat loads and split fractions varied in a single step; Approach 3: Two-level approach: network pinch varying both heat loads and split fractions

Table 5

 5 Split fraction of branches for pinched HEN by the new two-level approach (with respect to main streams)

	Branch No	15 (Str 10)	18 (Str 1)	21 (Str 9)	24 (Str 2)	27 (Str 3)
	Split fraction	0.472	0.850	0.682	0.534	0.452

Table 6

 6 Redistributed heat loads for pinched HEN by the new two-level approach

	Hx no	Redistributed heat loads (MW)	Existing heat loads (MW)
	1	11.87	13.25
	2	14.12	13.19
	3	0.93	0.86


Bálint Sinkovics
email: sinkovics@ett.bme.hu

Péter Gordon

Gábor Harsányi

Computer modelling of the laser ablation of polymers

Keywords: ablation of polyimide, laser ablation, thermal modelling

Laser ablation of polymers has been the subject of extensive studies for nearly three decades. Results of these investigations have shown that thermal effects play an important role during frequency tripled nanosecond Nd:YAG laser ablation of polyimide. Despite of the studies, effect of these processes has not been satisfactorily clarified. Thermal effects, photochemical and other phenomena still exist and so the thermal modelling of laser ablation is a very specialised problem. The investigation in this paper deals with special conditions in thermal models, which highly differ from more common conditions. Although there are numerous excellent thermal modelling software available, mainly in the field of electronics technology, these are not suitable or ready for handling special conditions, for example the removing of material from the domain under investigation, which is one of the most important areas. Hence we decided to develop a special thermal modelling tool. The model designing aspects which define the development process are presented in this paper.

Introduction

Despite of the wide use of laser technologies, interaction between the laser beam and materials has not been clarified in every detail so far, especially in the case of polymers [START_REF] Lippert | Chemical and Spectroscopic Aspects of Polymer Ablation: Special Features and Novel Directions[END_REF], [START_REF] Dyer | Excimer laser polymer ablation: twenty years on[END_REF]. Research results have shown that laser-material interaction is photothermal in the IR and photochemical in the deep UV region. In the case of polyimide thermal effects are pronounced in the ablation process [START_REF] Lippert | Fundamentals and applications of polymers designed for laser ablation[END_REF]. However, experiments in this field have shown that thermal and chemical effects appear at the same time, and many other physical phenomena also play important role.

In electronics technology the use of laser processing is very widespread and is gaining more and more importance due to the miniaturisation. In manufacturing, applied laser parameters (pulse energy, pulse repetition frequency, speed of the moving laser spot) are adjusted to the proper value usually heuristically depending on the work in progress. The main reason for using this approach is in connection with thermal effects appearing during the laser processing. Since etch rate of the laser depends on temperature, the material structure of the sample and accumulation of heat influences etch rate. In consequence, to control the laser processing, thermal distribution of the sample has to be determined, which can be made by thermal modelling of the laser ablation. By using such modelling tool, proper laser parameters can be determined easier and faster.

Nowadays, several excellent and usable 2D and 3D thermal simulation tools are available, but usually these are specialized in general thermal problems. Capabilities of the most important types of software tools are also examined, but thermal modelling of the laser process is a very special application. Among others geometry of the domain under investigation varies during the simulation, because laser shots remove material from the sample. This one condition alone precludes the possibility of applying general thermal simulation software.

The reasons mentioned above led to the development of a special thermal modelling tool. Thermal models for modelling laser ablation have already been presented in our previous works on the subject [START_REF] Gordon | Sinkovics: Thermal Simulation of UV Laser Ablation of Polyimide[END_REF]. However, in this paper we present a considerably improved modelling tool which is capable of handling nonuniform mesh and capable of applying real, measured laser beam distribution to each of the simulations.

Our software has been developed by Matlab. This solution is not optimal from the viewpoint of time consumption of the simulation, because there are faster algorithms exist, e.g. the Sincos transformation and many others [START_REF] Székely | THERMAN: a thermal simulation tool for IC chips, microstructures and PW boards[END_REF], [START_REF] Rencz | A fast algorithm for the layout based electro-thermal simulation[END_REF], [START_REF] Song | Numerical simulation of thermal loading produced by shaped high power laser onto engine parts[END_REF]. However time consumption of simulations was acceptable for our purposes.

Basic theory for modelling of laser ablation

Final goal of our work is to obtain a model for the laser ablation which can be built into the control software of a laser micromachining system. For this reason, our model approaches the laser-material interaction from a phenomenological viewpoint and does not deal with molecular details of the process.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

2 Generally, there are two types of approaches to laser processing models in relevant literature, one of them deals with atomic level of the processes. Aim of our models is to clarify the role of different physical phenomena in the laser-material interaction. These models also study, among others, details of the photothermal and photo-chemical effects [START_REF] Illy | Enhanced polymer ablation rates using highrepetition-rate ultraviolet lasers[END_REF], [START_REF] Yingling | Coarse-Grained Model of the Interaction of Light with Polymeric Material: Onset of Ablation[END_REF], [10] [11]; plume; plasma formation and their shielding influence [START_REF] Misra | Diagnostics of laser ablated plasmas using fast photography[END_REF], [START_REF] Bogaerts | Laser ablation for analytical sampling: what can we learn from modelling[END_REF], [START_REF] Zeifman | Multiscale simulation of laser ablation of organic solids: evolution of the plume[END_REF], [START_REF] Itina | Mechanisms of small clusters production by short and ultra-short laser ablation[END_REF]; appearing mechanical stresses and acoustic waves due to the explosion-like processes during the laser shot [START_REF] Allmen | Laser-beam interactions with materials 2nd[END_REF]; or, coupled phenomenon like the Marangoni-convection [START_REF] Lei | Numerical analysis of the competetive influence of Marangoni flow and evaporation on heat surface temperature and molten pool shape in laser surface remelting[END_REF]. Since these models are near to the laser-material interface, they are very complex and do not deal with the influence of the overall structure of the sample, e.g. how the sample directs the heat away from the placement of the laser shot. Our other group of models deals with the overall effect of laser processing, and ignores nearly all micro-level and other factors, e.g. the shape of the laser beam or structure of the sample when made from inhomogeneous material [START_REF] Yung | 355 nm Nd:YAG laser ablation of polyimide and its thermal effect[END_REF], [START_REF] Yung | XPS investigation of Upilex-S polyimide ablated by 355 nm Nd:YAG laser irradiation[END_REF], [START_REF] Yung | High repetition rate effect on the chemical characteristics and composition of Upilex-S polyimide ablated by a UV Nd:YAG laser[END_REF], [START_REF] Wu | Shin: A simplified predictive model for high-fluence ultra-short pulsed laser ablation of semiconductors and dielectrics[END_REF], [START_REF] Korfiatis | Numerical modelling of ultrashort-pulse laser ablation of silicon[END_REF], [START_REF] Chakraborty | The effects of turbulence on molten pool transport during melting and solidification processes in continuous conduction mode laser welding of copper-nickel dissimilar couple[END_REF].

Our model is a combination of models mentioned above due to the goals we have set. To consider the inhomogenities in thermal structure of the sample and the influence of the cumulative heat, our solution deals with the thermal distribution in the whole sample, not only that near to the laser-material interface. Besides the phenomenological viewpoint, our modelling tool makes it possible to take into consideration micro-level effects, shielding of the plasma, proportion of photo-thermal and photo-chemical phenomena, or other effects.

Essence of our solution is that energy of the laser shot is transformed into heat in the model and all necessary calculations are performed with the aim of thermal distribution [START_REF] Gordon | Sinkovics: Thermal Simulation of UV Laser Ablation of Polyimide[END_REF]. As our model does not deal with the details of different phenomena arising during the very short time of the shot itself, the transformation of the energy is calculated with no time consumption.

Theoretical background of the thermal cell method [24]

Principle of the thermal cell method is that cells are defined in the geometry, and temperature of those cells scrutinized. Cells are considered as infinitesimal objects located in the geometric centre of elements with finite volume. Thermal resistance between the cells are provided by the geometry of finite-volume elements and the thermal parameters of materials. Those parameters also determine thermal capacitance assigned to each cell. The illustration of thermal cells can be seen in Figure 1. Thermal resistance between cells is calculated by the thermal Ohm's law (following equations represent onedimensional case):

A x R ⋅ = λ ( 1
)
where R is thermal resistance, x is the distance between two cells, λ is the coefficient of thermal conductivity and A is the surface, which heat flux flows through. If thermal resistance is given, then the heat flux flowing between two cells is:

R T T Q n n n 1 + - = (2)
where T n is the temperature of the n th cell.

In the case where values of thermal resistance and capacitance assigned to cells are given, temperature of a cell in function of time can be calculated by the Taylor expansion of ( )

dt t T n + : ( ) ( ) ( )( ) ... ! 2 ' ' ! 1 ' ) ( 2 + + + = + dt t T dt t T t T dt t T n n n n (3)
The first derivative of the temperature can be originated from the energy balance equation of the cell:

( ) dt Q C t T t n n n ∫ = 0 1 (4)
where C n is the thermal capacitance of the n th cell:

V c C n ⋅ ⋅ = ρ (5) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 3
where c is the specific heat capacity, ρ is the density and V is the volume of the element. In equation ( 3) the higher order terms can be neglected, since we hold the value of the time step small, so the temperature of a cell in function of time can be approximated as:

dt C Q t T dt t T n n n n + = + ) ( ) ( (6)
Heat loss has to be considered in the cells on the boundary of the range. In heat flux calculations the Newton formula [START_REF] Song | Numerical simulation of thermal loading produced by shaped high power laser onto engine parts[END_REF] or the Stefan-Boltzmann law (8) apply:

) ( n a c T T A h Q - ⋅ ⋅ = (7) 4 0 n r T A Q ⋅ ⋅ ⋅ = σ ε ( 8
)
where h is the heat transfer coefficient, ε is the emissivity, 0 σ is the black body constant and a T is the temperature of the gas or liquid.

Description of light energy transformation to heat

Energy density and fluence distribution of a Gaussian laser beam can be written as follows:

        ⋅ + + - ⋅ = z y x e F z y x F α σ 2 2 2 2 0 ) , , ( (9)
where α is the absorption coefficient, x, y, z are the space coordinates. This means that the origin is at the centre of the beam and axis z is perpendicular to the surface of the material. The first part of the exponent describes the Gaussian distribution of the laser beam and the second comes from the Lambert-Beer Law of absorption. Parameter of width (σ) comes from the "minimal spot diameter" equation:

D fM d γ σ 2 min 44 . 2 4 = = (10)
where d min is defined as the position where the intensity drops to 1/e 2 , f is the focal distance, D is the beam diameter, M 2 is a beam quality parameter which is 1 for the TEM 00 Gaussian beam, and γ is the wavelength [START_REF] William | Laser material processing[END_REF].

To include all micro-level effects, modelling tool applies factors that determine the proportion of the beam energy transformed into heat. These factors have two well-separated roles. All effects that decrease the energy of the beam before it is absorbed by the material (i.e. plume formation, reflection, etc.) are considered in the coupling efficiency factor (CEF). Therefore, CEF is the ratio of energy coupled in the material and energy of the beam. Then energy is absorbed exponentially by the material, and transformed into heat by an efficiency factor we refer to as transform efficiency factor (TEF). TEF gives good estimation on the ratio of thermal and chemical processes that take place during the ablation process. Detailed description of how these factors are determined from real experiments was also presented in our previous work [START_REF] Gordon | Sinkovics: Thermal Simulation of UV Laser Ablation of Polyimide[END_REF].

F 0 is the fluence at the centre of the beam and it can be calculated from the pulse energy

) ( imp E : 2 0 2πσ imp E CEF F ⋅ = (11)
since the pulse energy is the definite integral of the fluence over the whole surface and it has to be decreased by CEF:

2 0 2 0 2 2 2 2 πσ σ F dxdy e F E CEF y x imp = = ⋅ + - ∞ ∞ - ∞ ∞ -∫ ∫ (12)
Energy (E) in a finite element of volume (its volume is

n m l × ×
) can be calculated based on (11):

n m l c E TEF T ⋅ ⋅ ⋅ ⋅ ⋅ = ∆ ρ ( 14
)
where c is the specific heat capacity, ρ is the density and TEF as described previously.

M A N U S C R I P T A C C E P T E D

Background of the thermal modelling tool

To summarize, the main object of the Matlab-based modelling tool is to solve equation ( 5) considering the geometry and material parameters, the initial and boundary conditions.

The flowchart of the solution can be seen in Figure 2. At first, thermal cells and connection between cells are created according to the geometry of the domain under investigation, the output is a special descriptive structure. This step is detailed in the next section. The second step is to increase the mesh density. The third step is to create a one-dimensional temperature vector, which gives the temperature of the thermal cells. At this stage of the solution this vector is filled up with the initial temperature of the material. The next step is to calculate the effect of the laser shot according to the principles presented in section 3, this step is showing the new temperature distribution. With this knowledge and knowing the given threshold temperature, portion of material that is ablated can be determined. In the next stage, the descriptive structure has to be refreshed. Cells correspond to a piece of material which is ablated and are removed from the descriptive structure. This means that the neighbouring relations between the thermal cells have to be updated, too. In the next stage of the process, the temperature vector is re-generated according to the new descriptive structure, since the number of thermal cells is decreased due to the ablation. Thermal simulation, which determines the thermal behaviour of the domain under investigation, takes place as last step. In the case of consecutive shots, the result of the simulation refers back to the fourth step. Outputs of the whole process are temperature distribution and the surface profile.

Definition of cells and thermal connections between cells

Definition of cells and calculation of the assigned geometric parameters are completed according to the available geometric data. The domain under investigation is divided into rectangular solids and at the beginning of the process its density is constant in each direction.

Calculation, storage, determination of placing cells within the grid and relation between neighbouring cells are also important points. A descriptive structure is created by the software to handle the cells and their relationships between neighbouring cells. This is shown in expression [START_REF] Itina | Mechanisms of small clusters production by short and ultra-short laser ablation[END_REF].

                          - + - + - + - + - + - + - - + + - - + + - - + + L ] [ ] [ ] [ ] [ 2 2 2 2 2 2 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 2 2 1 1 1 1 1 z z y y x x z z y y x x z z z z y y y y x x x x h h h h h h h h h h h h N N N N N N N N N N N N z y x C z y x C λ λ (15)
This structure is a special two-dimensional matrix. Data in columns are assigned to the elements of a thermal cell with a given ordinal number as a row vector, so that the elements of the structure are vectors. Elements of the first row are material parameters and therefore geometric data can be assigned to the given thermal cell. These data are sequentially: C -thermal capacitance, λ -coefficient of heat conductivity, x, y and zgeometric sizes of the element. The following six rows contain the ordinals of the neighbouring cells divided into rows per direction.

If the ordinal of the n th cell is included in one of the rows in the m th column (in vector N), it means that the n th cell is the neighbour of the m th cell from a specified direction, this direction defined by position of the ordinal in vector N. The vector in the last row contains the coefficient of heat transfer in six directions belonging to the specified element.

Creating a system of solvable equations

If descriptive structure is made as above, a system of differential equation has to be created, according to equation [START_REF] Székely | THERMAN: a thermal simulation tool for IC chips, microstructures and PW boards[END_REF], where the actual cell is denoted by m and its neighbours are denoted by n:

∑ ∑ ∑ ∑ - + - + - + - + - + - + ⋅ ⋅ - ⋅ + ⋅ + ⋅ - = z z y y x x N m n m m amb n m n z z y y x x N n m n m m n m m m dt R C t T t T h dt R R C t T t T dt t dT , , , , , , ) ( ) ( ) ( ) ( ) ( ) ( (16)
where T m (t) is the temperature of the m th thermal cell at moment t, C m is the thermal capacity of the m th cell, T n (t) is the temperature of the neighbours of the m th thermal cell selected according to the vector N, m n R is the thermal resistance of the m th cell in the direction defined by the vectors of N, n R is the thermal resistance of the neighbouring element defined by the vector N, m n h is the heat transfer coefficient in the specified direction, T n amb (t) is the ambient temperature in the direction specified. First term of the expression means that whilst the elements of the six vectors of contiguity are swept, the neighbours of the m th element in the direction specified are selected and equation ( 2) is applied to them. The second part of the formula is applied to the elements that have no neighbours, since they are located at the boundary of the region or are in connection with an ambient atmosphere. For these cases Newton formula is used. Elements with no neighbours are denoted by m N . Let us define multiplicators of temperature-change as follows. For the first part of the formula:

dt R R C e n m n m ⋅ + ⋅ = ) ( 1 1 (17)
For the second part of the formula:

dt R C h e m n m m n ⋅ ⋅ = 2 (18)
The definition of multiplicators in the case of neighbours in direction x is the following (direction y and z is similar):

dt x x C z y z y e m m n n m m m n n ⋅         ⋅ + ⋅ ⋅ ⋅ ⋅ = λ λ 2 2 )) ( ), min(( 1 (19)
The numerator contains calculation of common surface of two neighbouring elements. This formula defines the surface of common side in case of both elements. The smaller value is considered by using min() function. It enables to define common surface of elements, even if they have only one neighbour in the specified direction. Denominator contains thermal resistance, thermal capacitance of the specified thermal cell and thermal resistance of its neighbouring element. Multiplicator regarding heat transfer is calculated as follows:

dt x C z y h e m m m m m m n ⋅ ⋅ ⋅ ⋅ ⋅ = λ 2 ) ( 2 (20)

Refreshing descriptive structure

After a laser shot, a new temperature distribution is calculated according to the cumulative heat and the energy of the shot. By knowing the new distribution, the material parts that have higher temperature, than the threshold, are ablated. In this way a new geometry is evolved. Therefore, descriptive structure has to be refreshed. Cells which have higher temperature, than the threshold, are removed from descriptive structure. For this reason, the neighbouring conditions where cells are removed have to be modified: these element no longer have neighbours in the given direction, they are instead in connection with the atmosphere.

Modelling of laser ablation

Modelling parameters were consistent with the parameters of the Coherent AVIA 355-4500 frequency tripled, Q-switched Nd:YAG laser system, which is the appliance used in our laser experiments. Width parameter was 6 µm, according to the laser system. Material parameters correspond to parameters of 1 , which is a widely used polymer in electronics technology and serves as the basis of our laserpolymer investigations. Thermal parameters were the following [START_REF]MatWeb material database[END_REF]:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6 polyimide
K m W ⋅ = 29 . 0 λ , 3 1470 m kg = ρ , K kg J c ⋅ = 1130 (21)
According to our previous works, 0.132 TEF*CEF factor and 535 ˚C ablation threshold temperature and 1.4 m µ 1 absorption coefficient were chosen for performing the simulation [START_REF] Gordon | Sinkovics: Thermal Simulation of UV Laser Ablation of Polyimide[END_REF]. The initial temperature of the material was 25 ˚C. Initial mesh structure can be seen in Figure 3, size of mesh elements near the laser spot was 62.5 x 62.5 x 31.25 nm. For solving the time-dependent equations built-in Matlab functions (ode()) were used. These solvers apply adaptive time stepping, in our simulations the order of the minimal time step was 0.1 ns. Due to efficient solver algorithms, time consumption of the simulations was of the order of minutes2 per laser shots.

The result of simulation of a single shot with 10 µJ pulse energy is presented in Figure 4, the temperature of thermal cells is shown. In Figure 4 a) thermal distribution can be observed immediately after the laser shot and the ablation. Gaussian distribution of the laser shot and gradual absorption of the energy in layers of thermal cells can be seen. Besides thermal distribution, on top of the domain under investigation the ablated hole can also be observed. In Figure 4 b) temperature of cells is presented 10 ms after the laser shot. It can be observed that thermal energy of the laser shot has spread within the material, and the material cools down due to the heat loss in the domain under investigation.

Besides temperature distribution of the domain, other output of the modelling tool is the calculated surface profile. Since the tool is capable of handling any fluence distribution, measured distribution was applied in the simulation to verify the calculated profile. Fluence distribution of the laser system was measured by scanning with a high-speed silicon photodiode3 masked by an 5 µm pinhole 4 , the fluence distribution can be seen in Figure 5. In case of investigation of one laser shot, AFM5 micrograph of the ablation hole constituted the basis of verification for the simulation. Measured surface profile of the polyimide can be seen in Figure 6, where a single 18 µJ laser pulse was applied. Except fluence distribution, which is consistent with the measured values, applied modelling parameters are the same as in the previous case. The result of the simulation can be seen in Figure 7, where calculated surface profile in the symmetry plane of the laser beam is presented in comparison with the cross-section of the AFM micrograph. It can be seen that the simulated profile is in close proximity to the measured profile.

Consecutive shots are also simulated. Since perpendicular span of the AFM is limited, cross-section images of the polyimide foil constituted the basis of verification. In Figure 8 calculated surface profiles formed by punching can be seen. In the course of these experiments pulse repetition frequency of the laser was 30 kHz, and pulse energy was 22 µJ. By applying these parameters, considerable inhomogeneity appears in the beam, due to the frequency tripling of the laser system. Inhomogeneity can be observed in crosssection images and in simulation results as well, since the modelling tool is capable of handling the measured fluence distribution. As can be seen, simulated profiles approximate well to measured results, but it has to be noted that the modelling tool underestimates the ablated material in these experiments. The elimination of this problem requires further investigations.

Conclusion

A Matlab-based tool for the modelling of laser ablations of polymers from a thermal viewpoint is presented in this paper. Since the latter is a very special problem, a special approach has been applied to solving the equation of heat conduction. The finite difference method has formed the basis of the solution, a simple and expressive method. However, a novel descriptive structure has also been devised, which describes the thermal connections between the elements of the applied mesh. This descriptive structure has the following properties and abilities:

M A N U S C R I P T A C C E P T E D

1. Mathematical background of the solution is evident, and it can be implemented in Matlab or in other high-level programming language within a short time frame. 2. In Matlab environment, solver algorithms are already available and do not have to be implemented. 3. Due to the composition of descriptive structure, the modelling tool is capable of handling nonuniform mesh. With the use of such non-uniform mesh, accuracy of calculations is much more precise, high thermal gradient near to the laser spot can be handled easier and the time consumption of the simulation process decrease. 4. The descriptive structure makes it possible to apply further physical phenomena to take into consideration micro-level effects, such as plasma shielding or other side-effects. 5. The descriptive structure also makes it possible to handle any fluence distribution, so that real, measured distribution of the laser system can be applied in simulations. 6. The descriptive structure makes it possible to vary the geometry of the domain under investigation during the simulation. Comparison of simulated and measured profiles of a single and multiple laser shot have been presented as verification of our method. The next phase of our work will be directed at further investigations into the application of the modelling tool.

  thus the rise of temperature in a finite element of volume n m l × × can be calculated as:

Figure 1 :

 1 Figure 1: Thermal cells Figure 1 a): Connection of the thermal cells in the domain Figure 1 b): Nomenclature of thermal resistances

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 23456 Figure 2: A flow chart showing the most important steps in the Matlab-based modelling tool Figure 3: Mesh structure of the domain under investigation Figure 3 a): x-y plane, the placement of the laser spot is at the centre of this plane Figure 3 b): x-z plane, the placement of the laser beam is in the middle-top of this plane Figure 4: Result of simulation, temperatures of the thermal cells Figure 4 a): Immediately after the laser shot and the ablation Figure 4 b): After 10 ms following the laser shot Figure 5: Fluence distribution of Coherent AVIA 355-4500 frequency tripled, Q-switched Nd:YAG laser system (pulse energy: 18 µJ) Figure 6: AFM micrograph of ablation holes in Upilex-S, etched with one pulse of 18 µJ

Figure 7 :Figure 8 :

 78 Figure 7: Comparison of measured and simulated surface profiles (in the symmetry plane of the laser beam) Figure 8: Cross-section images of polyimide foil and simulated surface profiles

Figure 8 a

 8 Figure 8 a): after 5 shots Figure 8 b): after 10 shots Figure 8 c): after 15 shots Figure 8 d): after 20 shots

Upilex-S, high temperature polyimide film, product of UBE Industries Ltd.

Platform: Intel® Core™ 2 Duo 2 GHz with 2 GB of RAM

Purchased from the Thorlabs GmbH, part number: FDS010

Purchased from the Thorlabs GmbH, part number: P5S

Veeco Innova surface probe microscope was applied in contact mode