M Niklas

M Favre-Marinet

D Asendrych

Numerical simulation of microchannel network with complex geometry

Keywords: integrated circuits, cooling systems, microchannels, flow resistance

The paper presents the results of a numerical study devoted to the hydraulic properties of a network of parallel triangular microchannels (hydraulic diameter D h = 110 µm). Previous experimental investigations had revealed that pressure drop through the microchannels system dramatically increases for the Reynolds number exceeding value of 10. The disagreement of the experimental findings with the estimations of flow resistance based on the assumption of fully developed flow were suspected to result from the so-called scale effect. Numerical simulations were performed by using the classical system of flow equations (continuity and Navier-Stokes equations) in order to explain the observed discrepancies. The calculations showed a very good agreement with the experimental results proving that there is no scale effect for the microchannels considered, i.e. the relevance of the constitutive flow model applied was confirmed. It was also clearly indicated that the excessive pressure losses in the high Reynolds number range are due to the secondary flows and separations appearing in several regions of the microchannel system.

Introduction

The fast development of high-density (VLSI) microelectronic circuits requires more and more effective ways of cooling the microchips. In principle the higher the level of integrity (responsible for the CPU speed) the more heat is released as shown schematically in Fig. 1 prepared on the basis of data taken from [START_REF] Mahajan | Emerging directions for packaging technologies[END_REF]. The tendency presented here makes it clear that further progress in VLSI systems development is limited by the capabilities in heat removal from microprocessors. As a possible solution microchannel-based circuits with liquid cooling media can be employed as providing very high heat-transfer coefficient. Although such systems are being already constructed (see sample design in Fig. 2) their performance is still far from being quite satisfactory and a lot of research efforts have to be done to solve the existing problems.

The main advantage of small heat exchangers is their high cooling performance. Dongqing Li [2] showed that it may be even up to 60 times higher than the corresponding values for conventional, macroscale heat exchangers. The application of microchannels leads, however, to significantly increased flow resistance which is a serious drawback of this technology. Several early experimental studies provided results on the friction factor which cannot be easily explained with the use of conventional transport theories. Wu and Little [START_REF] Wu | Measurement of friction factors for the flow of gases in very fine channels used for microminiature joule-thomson refrigerators[END_REF][START_REF] Wu | Measurement of the heat transfer characteristics of gas flow in fine channel heat exchangers used for microminiature refrigerators[END_REF] investigated trapezoidal microchannels (etched in silicon and glass) and found the friction factor up to 30% higher than the value predicted for a circular pipe on the basis of conventional law. Similarly the work of Qu et al. [START_REF] Qu | Pressure-driven water flows in trapezoidal silicon microchannels[END_REF] revealed a significant increase of flow resistance in the wide range of Reynolds number studied. On the other hand recent data available (e.g. work of Baviere et al. [START_REF] Baviere | Experimental of water flow through smooth rectangular microchannels[END_REF]) do not reveal any deviations from the theoretical estimations even for very small-scale microchannels. An important aspect of microfluidics is, as reported in several papers, the premature (in terms of Reynolds number) transition to turbulence. As indicated by Peng et al. [START_REF] Peng | Frictional flow characteristics of water flowing through rectangular microchannels[END_REF][START_REF] Peng | Flow boiling through V-Shape microchannels[END_REF] the Fanning friction factor exhibited high levels for Reynolds numbers within the range of 200÷700 (depending on the hydraulic diameter of the microchannel). The authors interpreted this unexpected increase in f as the result of early transition.

As can be seen from the above review of the literature data, the knowledge and understanding of the fluid flow and heat transfer behaviour in microchannels are far from being complete. Some authors try to explain the discrepancies between experimental and theoretical results with the so-called scale effect suggesting the irrelevance of the conventional model of continuum for very small scales. There is a common feeling of the need for further fundamental investigations in order to provide more precise experimental data (due to small scales the measurement uncertainty is much higher than for macroflows) and look for the possible improvements of the constitutive models describing the microflows.

The present numerical study was undertaken within the framework of a larger research project devoted to the new design of a micro heat pipe with a system of triangular microchannels as a heat exchange section [START_REF] Mazellier | Etude experimentale d'ecoulements liquides en microcanaux[END_REF]. The experimental part of the project, preceding the current simulation work, allowed to find that the flow resistance through the network is unexpectedly high although still in the laminar flow regime [START_REF] Mazellier | Etude experimentale d'ecoulements liquides en microcanaux[END_REF]. In order to provide the complementary knowledge to the experimental trials suggesting the existence of a scale effect, the numerical simulation exactly followed the geometry of existing facility. It was believed that using conventional laws of hydrodynamics it would be possible to discover the sources of increased pressure losses and either to negate the scale effect or to support its existence.

Experimental arrangements

2.1. Network geometry and instrumentation. The test section under consideration consisted of a system of 55 parallel microchannels (see schematic view of the network in Fig. 3a) which were etched in a silicon wafer by a wet anisotropic technique and then covered by a silicon plate. The microchannels system was supplied with the circulating working liquid (demineralised water) by the external closed-loop circuit. The fluid was distributed among the microchannels through the inlet channel termed as distributor and then collected by the collector. The microchannels were of triangular cross section (see Fig. 3b) characterised by the hydraulic diameter D h = 110 µm and the length L = 18.4 mm so it corresponded to very long channels (L/D h = 167). The x axis of the coordinating system was placed along the symmetry axis of the system with its origin located in the front wall of the distributor (see Fig. 3a). The entire dimensions of the microchannels network were determined by its longitudinal (along x axis) and transverse (along z axis) extents h and w, respectively.

The pressure drop along the system ∆p being the main parameter of interest was determined by the difference of pressure at the inlet p in and outlet p out (∆p = p in -p out ) which were measured by means of two piezo resistive strain gauge transmitters placed upstream/downstream of the tested cell. The system bulk velocity was determined from the mass flow rate which was measured by precision weighing or by a highaccuracy flowmeter. The Reynolds number based on the average bulk velocity across the network U net and the hydraulic diameter of a single channel D h (Re = U net • D h /ν) was varied in the range 1 ÷ 700. The pressure losses allowed to determine the Fanning friction factor defined as

f = ∆p 2ρU 2 net D h L (1)
and then the Poiseuille number

Po = f • Re (2)
being the main parameter characterising the hydraulic performance of the network. 2.2. Head losses measurements. Measurements of the pressure drop across the system of microchannels [START_REF] Mazellier | Etude experimentale d'ecoulements liquides en microcanaux[END_REF] represented by the Poiseuille number are reported in Fig. 4 as a function of Reynolds number and compared to the theoretical prediction. One can note a dramatic increase of the flow resistance above the level corresponding to fully developed conditions (Po = 13.3) for the Reynolds number exceeding about 10. One of the possible explanations of that growth of friction factor could be the laminar-turbulent transition taking place inside the microchannel as indicated by Peng et al. [START_REF] Peng | Frictional flow characteristics of water flowing through rectangular microchannels[END_REF] who found experimentally that it can occur for much lower, than one would expect, Reynolds numbers (even below 700). It is not, however, the case of reported experiment as the continuous increase in the Poiseuille number starts for much lower Reynolds number values and for Re ≈ 300 , its value is roughly three times higher than the theoretical estimation. In order to find out the sources of these enormous discrepancies, a numerical simulation was undertaken taking into account all the probable physical phenomena occurring in the laminar single-phase flow.

Numerical simulations

3.1. Simulation strategy. Because of the complex geometry of the microchannels network it was decided to carry out the simulations in steps starting from the simplest case of a single microchannel and then extending the geometry up to the entire system following the existing experimental facility. Such an approach allowed to identify easily the consecutive sources of head losses and to quantify their relative importance versus Reynolds number. The consecutive flow configurations could be grouped into two classes, i.e. computations for single channel and a network of channels. In particular the following cases were considered:

• single microchannel (SM):

-fully developed flow -case SM1, -flow with entrance effects, i.e. with developing boundary layers -case SM2, -flow with local losses occurring at the inlet to the channel from the distributor -case SM3,

• network of microchannels (NM):

flow through the simplified network -geometry did not include the inlet to the distributor and the outlet from the collector -case NM1, -entire microchannels network -case NM2.

Numerical method.

The CFD (computational fluid dynamics) simulations of the flow under considerations were performed by means of the commercial code Fluent 5.4. The continuity (3) and momentum-transfer (4) equations for laminar flow regime were employed as the governing system of flow equations.

div U = 0 (3)

ρU j ∂U i ∂x j = - ∂p ∂x i + ∂τ ij ∂x j + ρg i . ( 4
)
In order to achieve the desired accuracy of calculations the equations were discretized using second order schemes with the second order upwinding method. The segregated solver was applied with the SIMPLEC pressure correction method.

Grid generation was done with the use of Gambit 1.3. The solution was treated as converged when the level of residuals was lower than 10 -5 . This was typically required to run appr. one thousand iterations for single channel geometry and several thousand iterations for a channels network. The grid quality was carefully checked by performing grid independence tests of the solution in terms of the pressure drop across the system. The consecutive calculations were performed for 3•10 5 , 6•10 5 and 9 • 10 5 volume elements and the greatest differences found were lower than 3%.

Results of the simulation. Single channel flow. Fully developed flow -case SM1.

As the starting point the simulation was conducted for the fully developed flow inside a single microchannel which could be achieved by applying the so-called periodic boundary conditions at the inlet and outlet cross sections. With this assumption no specifications were made about the flow-field parameters and the flow through the microchannel was driven thanks to the prescribed mass flow rate equal to the 1/55 of a total flux. For this simplest flow configuration the pressure drop computed from the recovered flow was only due to the friction at the bounding walls and no other effects influenced the hydraulic performance. The solution, as it could be expected, was Reynolds number independent and agreed excellently with the literature data (Bejan [START_REF] Bejan | Convection Heat Transfer[END_REF]), i.e. Po = 13.3. Flow with developing boundary layers -case SM2. The second step of numerical simulation was aimed at the analysis of the so-called entrance effect resulting from the increased flow resistance due to developing boundary layers downstream the channel inlet (see Fig. 5). In order to approach that flow configuration the "velocity inlet" boundary condition was employed at the inlet cross section with the uniformly distributed velocity deduced from the mass flow rate corresponding to the single microchannel. At the outlet from the computational domain the "pressure outlet" boundary condition was set where all the flow-field parameters (except for the static pressure assumed here as a reference value) are extrapolated from the interior. The results of the simulation are shown in Fig. 6, where they are put together with the experimental as well as the theoretically predicted value for reference case SM1, i.e. for fully developed flow. Noticeable growth of Po number appears for Re number greater than 300, while in experiment such a growth takes place already at Re = 10 and is much more intense.

One can conclude that the contribution of entrance effect into total head losses is not the crucial one but undoubtedly it should be taken into consideration during modelling of the flow.

Flow with local pressure losses -case SM3. The last flow configuration concerning single microchannel was targeted at estimation of local pressure losses appearing typically at the inlet to the duct, i.e. losses resulting from change of the flow direction and abrupt decrease of the cross-sectional area. The geometry under consideration was formed by the extension of geometry for former cases (SM1 and SM2) by a part of distributing channel corresponding to single microchannel, i.e. its transverse extent (along "z" axis) was equal to the spacing between channels of an entire network (see Fig. 7). The adequate simulation of such a flow case would require the application of only one inlet to the computational domain (e.g. "inlet 1" on one lateral side in Fig. 7) and two outlets (on the opposite lateral side denoted "inlet 2" in Fig. 7 and the outlet from the microchannel). As far as the single microchannel is concerned (without knowledge about flow-field nonuniformity along the distributor and collector) the correct definition of boundary conditions was not possible, so it was decided to consider the simplified model. The assumption was made about the symmetrical geometry, i.e. with application of two inlets ("inlet 1" and "inlet 2" in Fig. 7) where "velocity inlet" boundary conditions were used (similarly like in case SM2) with uniform velocity distribution deduced from the mass flow rate. At the outlet from the microchannel "pressure outlet" boundary condition was employed. As it can be found in Fig. 8 the Poiseuille number for case SM3 (open circles in Fig. 8) slightly exceeds data for flow with developing boundary layers (case SM2 -open triangles in Fig. 8) in the range of high Reynolds numbers (Re > 1000). The detailed analysis of the flow-field revealed the occurrence of separation regions just downstream the inlet to the microchannel and they are identified as being responsible for the Po number increase. The comparison of the numerical results with the experimental data presented in Fig. 8 clearly reveals that the entrance effect and the local losses have minor impact on the flow resistance and cannot justify the discrepancy in the range of moderate to high Re numbers.

Flow through the system of channels. The significant differences between numerical and experimental data in terms of head losses indicates the need of the simulation for entire channels network. Moreover, they seem to question the existence of the scale effect suggested to be the source of increased flow resistance.

The real system consists of the following elements (see Fig. 3a):

-55 triangular microchannels, -two transverse channels (distribution and collector channels) connecting the microchannels, -inlet and outlet from the system. At this stage of CFD simulation two flow configurations of microchannels network were investigated:

-simplified system (with neglected inlet ant outlet parts of the network), -actual network system including the inlet and outlet sections.

It was decided to distinguish the above cases in order to have a better insight into the physics of the flow and to quantify the independent sources of head losses in the network.

Flow through the simplified network -case NM1. The simplified model of a network [START_REF] Niklas | Numerical modeling of flow in microchannel[END_REF] is presented in Fig. 9. The following assumptions were applied:

-it was accepted that the model consists of even number of 54 microchannels, -only half of geometry (consisting of 27 microchannels) was considered due to its symmetry, -the inlet and outlet sections were not included in the geometry, -one of the ends of distribution channel was used as an inlet to the system (see Fig. 11) with "velocity inlet" boundary condition where uniform velocity distribution was set at the level resulting from the mass flow rate through the half-system, -one of the ends of collector channel was used as an outlet from the system (see Fig. 11) with "pressure outlet" boundary condition applied. The results of the simulation of the flow through the simplified system of microchannels are presented in Fig. 10 together with all previously obtained results. It is easily seen that the flow through the network introduces additional head losses which cause that the results (closed circles in Fig.

are much closer to experimental data (closed diamonds in Fig. 10) than those corresponding to the flow through single channel cases (cases SM1, SM2 and SM3). For the present configuration significant growth of Po number appears at Re number greater than 50 but the Po=f(Re) curve is not so steep like in experiment. Nevertheless, on the basis of these results it is possible to have an insight into the flow-field and try to answer the question what is responsible for so high growth of flow resistance in real conditions. Figure 11 shows the pressure variation along the distribution and collector channels corresponding to Re = 195. As it can be easily noted the static pressure in distributor is nearly constant (except the minor drop in close vicinity of the inlet) while in the collector a significant pressure drop is observed reaching 1/3 of the level corresponding to the farthest channel (z/w ≈ 0.5). The pressure field distribution can be explained with the use of the data of Fig. 12 where the velocity vector magnitudes are presented in both channels. While the velocity field in the distributor does not reveal any local nonuniformities (see Fig. 12a), the rapid "jets" outflowing from trian-gular microchannels are easily identified in the collector (see Fig. 12b). These jets resulting from an abrupt enlargement at the microchannels outlet have a great influence on the flow as they are responsible for intense mixing in the collector and for subsequent high head losses in this region. These additional head losses cannot be simply avoided because they are related to the sudden change of cross-sectional area (triangular microchannel -collector channel) and due to technological problems cannot be replaced by diverging channel endings. Flow through the microchannels system -case NM2. As the last step of the present numerical work the simulation was carried out for the computational domain exactly following the real microchannels network. Taking advantage of the symmetry of the system the computations were limited to the half-geometry [START_REF] Niklas | Numerical modeling of flow in microchannel[END_REF][START_REF] Niklas | Pressure losses in a network of triangular microchannels[END_REF][START_REF] Niklas | An experimental study and numerical modeling of the flow in a network of microchannels[END_REF]. The view of the geometry under interest including the inlet and outlet sections is shown in Fig. 13. The boundary conditions applied for the current configuration were analogous with those used for simplified network model (case NM2), i.e.:

-"velocity inlet" set at the inlet to the system where uniform velocity distribution was assumed, -"pressure outlet" imposed at the outlet from the network where no specifications about the flow-field were done. and additionally the "symmetry" condition was set at a plane of symmetry (see Fig. 13). The results of computations for entire network geometry shown in Fig. 14 are qualitatively identical with the experimental data revealing, however, the slight overestimation of Poiseuille number. The differences may be due to the uncertainties in the numerical simulation, but also to measuring errors, caused, for instance, by presence of air trapped in the corners of the triangular cross-section. According to the present numerical and experimental results, it seems that no size effect is present in the flow investigated with the characteristic lengthscale of 110 µm. It means that the constitutive model (Navier-Stokes) is able to adequately predict the flow in the network under consideration. Careful analyses of the flowfield in this complex system of channels show that there are some other sources of high pressure losses in addition to those typically appearing in microflows. Fig. 15 illustrates the pressure drop in the distributor and collector channels for the selected values of Reynolds number (Re = 3, 195, 390, 584) covering the experimental range. The results are presented here in normalised way in order to allow the comparison of pressure drops varying in a very wide range. The normalised pressure was defined as follows

p norm = p -p out p max -p out (5)
where p max represents the maximum pressure. As it can be found in Fig. 15 the highest pressure drop is observed at the inlet to the distribution channel (except for the case of the lowest Re) and moreover, this pressure drop is strongly Re number dependent (the higher the Re number the higher is the pressure drop). The explanation of such a pressure distribution is clearly visible in Fig. 16 presenting the contours of the z-velocity (transverse) component in the inlet region. The sudden change in the flow direction at the end of the inlet to the distributor gives rise to a strong separation just downstream the inlet corner. The coherent vortex attached to the distributor wall with extensive reverse flow (dark region in Fig. 16) is responsible for significant reduction of the distributor channel crosssectional area and acceleration of the stream leading in turn to the high pressure drop in this region. In the farther transverse distances exceeding the dimension of separation zone the pressure is practically constant. As it can be easily estimated from data collected in Fig. 15 the inlet separation region is responsible for more than 20% of total head losses and its contribution goes up for increasing Reynolds number. In order to interpret the pressure distribution in the collector presented in Fig. 15 the insight into a single microchannel can be of help. Fig. 17 shows the static pressure variation along the centreline of the system (z = 0). The initial part of the curve corresponds to the inlet section (x/h < 0) and to the distribution channel (as shown in the scheme at the top of Fig. 17). The constant pressure gradient observed in the channel is a convincing evidence of the fully developed flow. The nonlinear pressure distribution in the inlet part of the microchannel, caused by the entrance effect, is practically disappearing for x/h > 0.3. In the ending part of the p norm = f(x/h) distribution corresponding to the collector and the outlet section, a very high pressure drop is observed (x/h ≥ 1). The source of these losses is shown in Fig. 18 where the contour map of x-velocity (longitudinal) component in the outlet region is presented. Similarly like in the inlet region the sharp corner leading to sudden change of flow direction causes the origin of separation vortex (dark-coloured region in Fig. 18). It noticeably contributes to the total pressure drop and, as can be estimated from pressure distribution shown in Fig. 17, its contribution may also (like separation eddy in the inlet section) exceed 20%. Summing up both the separation zones in the microchannels network may be responsible for nearly half of the total head losses. An important aspect of the flow through the network of microchannels is the nonuniform distribution of the fluid flow rate among the channels. Fig. 19 shows the variation of channel bulk velocity U ch normalised by the network bulk velocity U net across the network for two selected values of Reynolds number, namely Re = 195 and 390. The mean velocity varies significantly from channel to channel as a result of the nonuniform pressure differences (see Fig. 15). The highest velocities and mass flow rates are obviously observed in the central part of the system with maximum/minimum values ratio about 3. It may be noted that the distribution of the bulk velocity along the system is slightly Reynolds number dependent (the higher Re number the higher (U ch /U net ) max value -compare the bulk velocity profiles in Fig. 19a and19b). It should also be remarked that such a flow rate nonuniformity across the network results in nonuniform heat removal capabilities and should be regarded as a drawback of a cooling system based on the present configuration. It is worth emphasizing that the objective of the research was to investigate the flow characteristics of a test section with the same geometry as that of a micro heatpipe. The test-section arrangement was obviously not optimized from the point of view of heat exchangers.

Summary and conclusions

As it was shown in the preceding sections of the paper the CFD simulation based on classical laws of hydraulics may be an adequate tool for flow analysis of a complex microchannels network. The computations clearly showed that there is no scale effect (or it is of negligible importance) in triangular microchannels with a hydraulic diameter of 110 µm. This research also underlines the care which must be taken in interpreting pressure measurements during the study of flow in complex microchannels systems. The flow resistance in one microchannels cannot be estimated on the basis of pressure drop in the entire network. It requires a careful separation of pressure variations in the region of interest and in the parts of the system which are used to connect the microchannel to an external hydraulic circuit. The importance of this separation is clearly visible in Fig. 20 where the contribution of the various phenomena (different sources of pressure losses in the flow) to the total pressure losses as a function of Reynolds number is shown. The curves correspond to the various assumptions (flow configurations) used for the computations starting from single microchannel and finishing on the entire microchannels system. The following phenomena were identified and quantified thanks to step-by-step simulation procedure:

-viscous losses -representing the losses due to friction observed in fully developed flow (case SM1), -entrance effect -being a measure of the pressure drop caused by developing boundary layers in the initial part of channels (case SM2), -local losses -generated in the inlet to the microchannels and associated to abrupt contraction (included in case SM3), -mixing -caused by jets exiting the triangular microchannels and entering the collector channel (analysed in case NM1), -recirculation -head losses appearing in the inlet and outlet sections due to separation of boundary layers as a result of improper geometry (case NM2).

As can be found from Fig. 20 the flow resistance is mainly dependent on three phenomena, i.e. viscous losses, mixing and recirculation, having jointly a major contribution to the total losses in the system analysed. The importance of friction losses is really huge (nearly 80%) for very low Re which prevents the generation of separation zones. For increasing Re numbers recirculation regions start to appear in the flow and viscous losses contribution is falling down to the level of appr. 20% for transitional flows although remaining at the same absolute level (Po = 13.3). The head losses due to separation zones as well as mixing in collector channel are Reynolds number dependent in similar manner and they are likely proportional to the flow rate squared. As it is indicated in Fig. 20 the total pressure drop in the range of high Reynolds numbers is dominated by losses in the distributor, collector and outlet channels. For Re higher than 200 the pressure drop caused by friction in the microchannels is lower than half of the total losses. It shows that optimization of the flow in microchannels systems should be concentrated on the proper design of the distribution and outlet channels (for the sake of separation) and of the collector channel (for the sake of mixing). The two remaining phenomena influencing the flow resistance, developing boundary layers and inlet losses are in general of minor importance as for Re = 700 (the highest Re number shown in Fig. 20) their common contribution is appr. 7% with 5%-contribution of entrance effect. Even for optimised network design (with significantly reduced mixing and separation losses) these sources of head losses would be much lower than the pressure drop due to viscous friction.

The current investigations showed that careful design of complex network of microchannels, as in micro heat exchangers, is necessary to avoid excessive pressure losses in the system. It was also shown that computational fluid dynamics can provide the engineers the powerful tool enabling for the optimisation of particular design before the real installation is constructed.

Fig. 1 .

 1 Fig. 1. Historical power trend for Intel CPUs (after Ref. 1)

Fig. 3 .

 3 Fig. 3. The geometry of microchannels system; entire view of the network (a); cross-sectional view of silicon wafer (b)

Fig. 4 .

 4 Fig. 4. Pressure losses in a system of microchannels -experiment versus fully developed flow (after Ref. 9)

Fig. 5 .

 5 Fig. 5. The flow pattern in a flow with developing boundary layers

Fig. 7 .Fig. 8 .

 78 Fig. 7. Velocity field for the case of single channel with local (inlet) losses

Fig. 9 .

 9 Fig. 9. The geometry of a simplified network

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Contribution of the mixing to the flow resistance of the microchannels system

Fig. 12 .

 12 Fig. 12. Velocity vector magnitude distribution; inside the distributor (a); inside the collector illustrating mixing losses (b)

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Half geometry of the actual network

Fig. 16 .

 16 Fig. 16. Z-velocity (spanwise) component distribution in the inlet zone

Fig. 17 .

 17 Fig. 17. Pressure variation along the centreline of the system

Fig. 19 .

 19 Fig. 19. Distribution of microchannel bulk velocity across the network for: Re = 195 (a); Re = 390 (b)

Fig. 20 .

 20 Fig. 20. Contribution of all pressure losses into total flow resistance in the network

Acknowledgements. The Region Rhône-Alpes is gratefully acknowledged for the financial support under the project "Micro heatpipes" carried out by a network of French research institutes. Additional funding was provided by Częstochowa University of Technology as Statutory Funds under the grant No BS-1-103-301/2004/P. The authors would also like to thank the European Commission for the support of Mariusz Niklas stay in Grenoble in frames of Socrates-Erasmus programme.