Nathan Richermoz
email: nathan.richermoz@cea.fr

Damien Gloriod
email: damien.gloriod@cea.fr

Thierry Baffie
email: thierry.baffie@cea.fr

Zoé Anxionnaz-Minvielle

Pierre Coste
email: pierre.coste@cea.fr.

Roughness effect on thermal-hydraulic performances of additively manufactured meandering mini-channels

, show that the Nusselt number is not a linear function of the friction factor.

 mentioned these specific AM different roughness features, such as balling, hatch spacing, keyhole, dross… In the case of mini-channels, the flow is even more sensitive to these features and to high relative roughness. To study the AM roughness, there are mainly two numerical approaches.

. This approach reproduces well the thermal-hydraulic performances of surfaces but needs an accurate surface characterization and a high simulation power. In fact, such surface meshing requires about 30 cells for each peak or between peaks. In other words, the flow around roughness peaks and through roughness valleys has to be captured by the grid. The technologies to characterize surfaces topologies are available at a laboratory scale but it requires open surfaces or else destructive measurements.

Introduction

The potential of Additive Manufacturing (AM) technologies to produce heat transfer devices has been recognized for at least fifteen years, with an application to novel heat sinks designs characterized by improved pressure loss and heat transfer [START_REF] Wong | Pressure Loss and Heat Transfer Through Heat Sinks Produced by Selective Laser Melting[END_REF], leading today to a wide range of materials and methods [START_REF] Nafis | Additive Manufacturing for Enhancing Thermal Dissipation in Heat Sink Implementation: A Review[END_REF]. Consistently with the idea of enhancing thermal dissipation in heat sinks, the idea of using AM technologies to novel heat exchangers designs has also received interest, taking advantage of the wide new manufacturing possibilities [START_REF] Koca | A Novel 3D Printed Air-Cooled Fuel Cooler Heat Exchanger for Aviation Industry[END_REF]. Nevertheless, this process is quite recent and the surface topologies are not perfectly known and controlled because there is a strong relationship between rough elements and printing parameters [START_REF] Thole | Review of advances in convective heat transfer developed through additive manufacturing[END_REF][START_REF] Wildgoose | Impact of Additive Manufacturing on Internal Cooling Channels With Varying Diameters and Build Directions[END_REF][START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF][START_REF] Kaur | State-of-the-art in heat exchanger additive manufacturing[END_REF]. Relative roughness values up to 38% of the mini-channel hydraulic diameter are reported in the literature [START_REF] Kirsch | Heat Transfer and Pressure Loss Measurements in Additively Manufactured Wavy Microchannels[END_REF].

Furthermore, even if many new geometries are feasible thanks to Laser-Powder Bed Fusion (L-PBF), the AM technology used in this article to build the mock-ups, one can observe a collapse phenomenon of the channel's ceiling [START_REF] Snyder | Build Direction Effects on Microchannel Tolerance and Surface Roughness[END_REF][START_REF] Khan | Parameters optimization for horizontally built circular profiles: Numerical and experimental investigation[END_REF][START_REF] Çalışkan | Investigation of direct metal laser sintering downskin parameters' sagging effect on microchannels[END_REF] or partial obstruction and resulting strong deviations of cross-section shapes, especially when considering mini-channels (up to one third deviation compared to the specified hydraulic diameter) [START_REF] Wildgoose | Impact of Additive Manufacturing on Internal Cooling Channels With Varying Diameters and Build Directions[END_REF], [START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF], [START_REF] Arie | Experimental characterization of an additively manufactured heat exchanger for dry cooling of power plants[END_REF][START_REF] Niknam | Additively manufactured heat exchangers: a review on opportunities and challenges[END_REF][START_REF] Da Silva | Thermal and hydrodynamic analysis of a compact heat exchanger produced by additive manufacturing[END_REF][START_REF] Dahmen | Characterization of channels made by laser powder bed fusion and binder jetting using X-ray CT and image analysis[END_REF][START_REF] Snyder | Understanding Laser Powder Bed Fusion Surface Roughness[END_REF]. As a consequence, discrepancies are observed between thermal-hydraulic performances of the CAD (Computer-Aided Design)-based component and the AM one as well as between identical prototypes made of different materials [START_REF] Arie | Experimental characterization of an additively manufactured heat exchanger for dry cooling of power plants[END_REF]. It may induce prediction difficulties. For instance, some authors show that the friction factor from their experimental data for straight mini-channel is more than two times higher than the one from Blasius law prediction [START_REF] Favero | Additive manufacturing for thermal management applications: from experimental results to numerical modeling[END_REF] or that it differs from Blasius or Colebrook laws according to the roughness-to-hydraulic diameter ratios [START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF]. As a consequence, for a given pressure ratio, the increase of friction factor may lead to a decrease (up to 70%, [START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF]) of the mass flow rate which might

The second approach is roughness modeling, where "modeling" has to be heard in the sense that roughness effects are averaged in sub-grid scale models. It means that in the computational domain, the surfaces are perfectly smooth in the meshing but the roughness is taken into account with special treatments of the walls. The work made by the author [START_REF] Aupoix | Roughness Corrections for the k-ω Shear Stress Transport Model: Status and Proposals[END_REF] stands in the modification of the turbulent Prandtl number boundary condition at the walls. The authors [START_REF] Mazzei | CFD Modelling Strategies for the Simulation of Roughness Effects on Friction and Heat Transfer in Additive Manufactured Components[END_REF] have then applied this model to the authors [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF] experiments. The authors [START_REF] Mazzei | CFD Modelling Strategies for the Simulation of Roughness Effects on Friction and Heat Transfer in Additive Manufactured Components[END_REF] propose two correlations to determine the constants needed in their model. Along the same lines, the authors [START_REF] Favero | Additive manufacturing for thermal management applications: from experimental results to numerical modeling[END_REF] modify the condition close to the wall by adjusting the Equivalent Sand Grain Roughness (ESGR), using the implemented models in Ansys Fluent v20.1, until an optimum is found. This study reveals a strong relationship between roughness and printing orientation. However, their geometries were too complex to deduce a correlation linking the ESGR and the surface orientation.

The present work deals with an AM application in the context of a heat exchanger/reactor intended for process intensification in the chemistry industry. Its design is based on wavy minichannels. This shape is chosen because of its high thermal performances with moderate pressure drop while triggering the laminar/turbulent transition at lower Reynolds numbers, which is of great interest for chemistry applications [START_REF] Anxionnaz | Heat exchanger/reactors (HEX reactors): Concepts, technologies: State-of-the-art[END_REF][START_REF] Anxionnaz-Minvielle | Influence of the meandering channel geometry on the thermo-hydraulic performances of an intensified heat exchanger/reactor[END_REF][START_REF] Anxionnaz | Transposition of an exothermic reaction from a batch reactor to an intensified continuous one[END_REF]. It enables to combine high residence time with heat and mass transfer intensification in a Reynolds number range between 200 and 6,000.

The objective of the present work is: in a first step, to analyze and compare experimentally different mock-ups obtained either by Hot Isostatic Pressing (HIP) or by L-PBF (AM), with different cross-sections (circular and square) and various surface states (raw and polished); in a second step, to propose a numerical approach of the roughness modeling kind and to validate it on present experimental data.

The paper is organized as follows. The experimental set-up is described. Then the experimental results are presented, considering thermal, hydraulic and global performances, comparing HIP versus AM, circular versus square cross-sections and as-built versus polished minichannels. The numerical study is then presented, starting with an emphasis on the geometry, ending with results and a sensitivity analysis.

Experimental set-up

Description of the facility

The schematic diagram of the experimental rig is presented in figure 1. Figures 2 and3 illustrate the test facility which is composed of two loops.

A hot water closed loop, which temperature is controlled by a thermostatic bath and which flow rate is controlled by a magnetic pump. The instrumentation includes two temperature sensors (Pt100) located at the inlet and outlet of the channel, a flow rate sensor and a pressure drop sensor between inlet and outlet of the mock-ups, as illustrated in figure 3.

A cold water open loop, which instrumentation includes: two temperature sensors (Pt100), for inlet and outlet of the cooling boxes and a flow rate sensor.

Description of the mock-ups

The mock-ups made of Alloy 600, were fabricated by L-PBF and tested in as-built state.

Although Alloy 600 has applications in chemical, heat treatment and aeronautical industries, few studies relative to L-PBF of this alloy are available in the literature [START_REF] Anderson | NEET-AMM Final Technical Report on Laser Direct Manufacturing (LDM) for Nuclear Power Components[END_REF]. Our team recently developed L-PBF process parameters of this alloy; we studied the microstructure and mechanical properties of as-built and heat-treated samples [START_REF] Salvador | Microstructure and Tensile Properties of Alloy 600 Parts Produced by Laser Powder Bed Fusion (L-PBF) Process[END_REF], [START_REF] Gloriod | Characterization of Alloy 600 heat exchanger-reactor mock-ups obtained by Laser Powder Bed Fusion (L-PBF) process[END_REF]. We showed that a 10-45 µm gas atomized Alloy 600 powder, with limited amounts of impurities such as boron, phosphorus and sulfur, was compatible with L-PBF process. Built samples were close to 100% density.

Several heat exchanger-reactor mock-ups were produced. Their external dimensions are 105 mm in length, 22 mm in height and 6 mm in thickness. All mock-ups include a 2D zigzag geometry mini-channel with a hydraulic diameter around 2 mm, with either square or circular crosssections. The hot water flows into this mini-channel and the cold water flows on the upper and lower faces of the mock-up. One pattern of the mini-channel geometry is depicted in figure 4. Its nominal dimensions are =135° (bend angle), L=3.5 mm (half straight length between tow bends) and R=1.5 mm (bend curvature radius). This pattern is repeated seven times to make a whole channel (figure 5).

As illustrated in figure 5, the channel pattern is drawn on the outer surface of each mockup so as to be located prior to mock-ups cutting. The inlet and outlet pipes have been welded after L-PBF step. The thickness of the mock-up's wall separating the hot and the cold fluid is 2 mm.

The cold fluid exchanges heat on the two mock-up's external surfaces of 13 mm width and 78 mm length and flows in an around 4 mm hydraulic diameter channel. Four different AM mock-ups are studied (table 1). A reference mock-up, obtained by HIP and characterized in a previous work ([44]), is introduced to compare the influence of the manufacturing process.

The circular cross-section mock-up is considered as the nominal one. The channel wall roughness was measured on a part of the wall which was vertically manufactured which corresponds to the smoother part. Depending on the measurement method (optical or mechanical) and on the location, some discrepancies are observed. Therefore a correction factor is applied to these values in our numerical model, as described in the roughness modeling section.

In this paper, the classical one-dimensional parameters Ra and Rz are used to characterize the surface based on the vertical deviations of the roughness profile from the mean line. Let zj be the altitude of a point part of a surface profile.

Ra is the arithmetic average roughness:

𝑅𝑎 = 1 𝑛 × |𝑧 -𝑧̅ | (1)
Rz is the maximum peak to valley height of the profile:

𝑅𝑧 = max 𝑧 -𝑧̅ (2)
Where

𝑧̅ = 1 𝑛 × 𝑧 (3)
Baseline 'circular" mock-up was cut longitudinally by electro-discharge machining for optical roughness measurement of the channel walls; measurement based on focus variation was done with an Infinite Focus microscope from Bruker-Alicona. Average roughness (Ra) and maximum roughness (Rz) were measured on vertical walls (when referring to the mock-up orientation on the L-PBF baseplate).

A second circular cross-section mock-up, 'Circular bis', is also presented. It has been obtained with different printing parameters as shown in table 1.

The second mock-up channel, 'Square', has a square cross-section (2 mm width). The purpose of this mock-up is to compare the thermal-hydraulic performances of the wavy minichannel obtained from AM, 'Square', with the ones from HIP, named 'Smooth square' in table 1.

The HIP process consists in plate laser machining. The plates are then stacked and assembled under high controlled temperature and pressure. Unlike AM, circular cross-section are not feasible by HIP. That is the reason why AM square cross-section mini-channel has been manufactured. The influence of mini-channel manufacturing process on thermal-hydraulic performances as well as the performances of circular cross-section vs square one are studied in the 'experimental results' section.

The last mock-up, 'Polished', has a circular cross-section like the first one, but a chemical surface treatment of the channel has been done (50°C, HNO3-HF, 1 hour). By this way, the roughness has been modified from Ra =7.7 μm to Ra=4.1 μm and from Rz=44 μm to Rz=27 μm. Once the channel was chemically milled, the mock-up was cut longitudinally. The roughness has been characterized on the vertical walls parallel to their direction using a direct roughness meter.

Several printing parameters have been used to solve various printing issues ( [START_REF] Gloriod | Characterization of Alloy 600 heat exchanger-reactor mock-ups obtained by Laser Powder Bed Fusion (L-PBF) process[END_REF]). Table 2 presents the printing parameters common to all mock-ups. Parameters that have been changed between the mock-ups are specified in table 1 (powder quality and baseplate temperature). The two powders used, virgin and re-used, come from the same batch. The re-used powder was employed three times in the L-PBF machine; it was sieved after each printing cycle to remove the coarse particles. We noticed that its flowability has deteriorated with the cycles, which may impact the powder bed density which, in turn, can generate defects such as pores in the final parts. Mock-ups are printed in such a way that the zigzag amplitude is in the vertical direction (Figure 5a).

Data reduction

The temperature difference between the hot fluid outlet temperature Toh and inlet one Tih is:

∆𝑇 = 𝑇 -𝑇 (4𝑎)
Similarly the temperature difference between the cold fluid outlet temperature Toc and inlet one Tic is:

∆𝑇 = 𝑇 -𝑇 (4𝑏)
Thus, the logarithmic mean temperature difference of the heat exchanger is calculated as:

𝐿𝑀𝑇𝐷 = (𝑇 -𝑇 ) -(𝑇 -𝑇 ) ln(𝑇 -𝑇 ) -ln(𝑇 -𝑇 ) (5)
The Reynolds number Re, the Prandtl number Pr and the Nusselt number Nu, are calculated as:

𝑅𝑒 = 𝜌 • 𝑣 • 𝐷 𝜇 (6) 𝑃𝑟 = 𝐶𝑝 • 𝜇 𝜆 (7) 𝑁𝑢 = ℎ • 𝐷 𝜆 ( 8
)
where 𝑣 (m•s -1 ) is the flow averaged velocity, Dh (m) the hydraulic diameter, h

(W•K -1 •m -2 )
the convective heat transfer coefficient of the cold or the hot fluid. The water material properties are calculated at the average temperature of the fluid: viscosity

μ (Pa•s), density ρ (kg•m -3 ), conductivity λ (W•m -1 •K -1 ), heat capacity Cp (J•kg -1 •K -1 ).
On the cold fluid side, the heat transfer coefficient is determined from the Colburn correlation:

𝑣 = 𝑄 𝐴 • 𝜌 (9) 𝐷 = 4 • 𝐴 𝑝 (10) 𝑁𝑢 = 0.023 • 𝑅𝑒 . • 𝑃𝑟 (11)
Where Q (kg•s -1 ) is the mass flow rate, A (m 2 ) the cross section area of the channel and 𝑝 (m) its perimeter.

The thermal power Pth (W) is written:

𝑃 = 𝑄 • 𝐶 • ∆𝑇 ( 12
)
Then the intensification factor UA/V (W•K -1 •m -3 ) of the heat exchanger is calculated as:

𝑈𝐴 𝑉 = 𝑃 𝑉 • 𝐿𝑀𝑇𝐷 (13)
V is the hot fluid volume contained into the heat exchanger mini-channel (m 3 ).

Then the heat transfer coefficient of the hot fluid hh (W•K -1 •m -2 ) is determined:

1 ℎ = 𝐴 𝑈𝐴 - 𝑒 • 𝐴 𝜆 • 𝐴 - 1 ℎ • 𝐴 𝐴 (14)
Where Ah (m²) is the area of the wetted surface of the hot fluid; Ac (m²) the area of the wetted surfaces of the cold fluid; Am (m²) the mean between these two last ones. e (m) is the minimum thickness between the two fluids.

𝜆 (W•K -1 •m -1 )
is the thermal conductivity of the heat exchanger material (Inconel 600) and hc

(W•K -1 •m -2 )
is the heat transfer coefficient of the cold fluid, determined by Eqs. ( 8) and [START_REF] Çalışkan | Investigation of direct metal laser sintering downskin parameters' sagging effect on microchannels[END_REF].

Finally, the Nusselt number of the hot fluid is calculated according to Eq. ( 8).

Considering the hydraulic analysis, the Fanning number is calculated as:

𝑓 = ∆𝑃 • 𝐷 2 • 𝜌 • 𝑣 • 𝐿 (15)
Where ΔP (Pa) is the pressure drop between the inlet and the outlet of the hot fluid and Ldev (m) is the developed length of the channel.

Sensor's uncertainties

Measurement uncertainties are summarized in table 3.

As presented in table 3, experimental uncertainties are limited, so that even if they are all included in the following graphs, they are not always visible because hidden by the point size.

Experimental results

Manufacturing process comparison: AM vs HIP

The square cross-section channels have been manufactured either by AM or HIP. Their thermal-hydraulic performances are compared in this section. Two different sets of experiments have been performed: adiabatic experiments for hydraulics performances and non-adiabatic ones for thermal performances. As shown in figure 6 and figure 7, the heat exchange through the channel increases when the flow rate increases as expected. The heat transfer is enhanced in the AM channel. At the same average flow velocity, the intensification factor of the AM mock-up is 30% higher than the one of the mock-up made by HIP. This trend is similar for the Nusselt number as a function of the Reynolds number, which allows an average increase of 80% with L-PBF compared to HIP. As reported in literature [START_REF] Chen | Numerical simulation of laminar heat transfer in microchannels with rough surfaces characterized by fractal Cantor structures[END_REF], [START_REF] Dharaiya | A numerical study on the effects of 2d structured sinusoidal elements on fluid flow and heat transfer at microscale[END_REF], the heat transfer enhancement measured in AM mini-channels is related to higher roughness values in the AM channel.

Thermal performances

Hydraulics performances

The hydraulic characterizations of mock-ups are displayed as a function of the flow velocity and Reynolds number (figure 8 As seen in figure 8 and in figure 9, both the pressure loss and the Fanning number are always higher (up to 37%) with AM than with HIP. This result is directly related to the roughness, which is higher in 'square' AM mock-up than in HIP 'square' one, as seen in table 1.

Global thermal-hydraulic comparisons

Finally, the thermal-hydraulic behavior of the mock-ups can also be analyzed considering a global criteria: the intensification factor, UA/V (W•K -1 •m -3 ) as a function of the volumetric pumping power (W•m -3 ) as illustrated in figure 10.

This graph can be seen as a picture of the "operating cost" (pumping power) to "return on investment" (heat exchange) ratio. The uncertainties reported in the figure are discussed in the following section. From this point of view, the AM manufacturing is in average 30% more competitive than the same geometry made by HIP. This conclusion is very positive for AM because it means that the process is of interest, not to mention the new design opportunities.

Different AM mock-ups comparisons

Cross-section scans

Computed tomography (CT) scans have been performed on some 'Circular' mock-up crosssections (figure 11). Analyses were performed under 220 kV using a helicoidal Nikon XTH225 device equipped with a 5 µm micro-focus X-ray source. Voxel sizes were in the 35-39 µm range.

These CT scans characterize the geometrical variation of the channel between CAD dimension (2mm in diameter) and the real printed channel. Results showed the effect of the part orientation on the walls quality; channel sections were slightly smaller than on CAD designs; ceiling surfaces were rougher and more defective than floor surfaces. Circular sections were ovalized; the greatest reductions in section were noticed at the bends built in cantilever position.

In an attempt to take it into account, a "geometrical factor uncertainty" of 10% is introduced on the hydraulic diameter. This "geometrical factor uncertainty" is called herein uncertainty because (i) the whole set of mock-ups has not been fully analyzed by a CT scan and (ii) different discrepancies from the CAD geometry could have been obtained with the different printing parameters. However, this factor can also be seen as a sensitivity to the hydraulic diameter on the results.

Thermal performances

The thermal performances of AM mock-ups are shown for the intensification factor (figure 12) and the Nusselt number (figure 13) as a function of the flow velocity and the Reynolds number respectively.

In figure 12 and figure 13Figure 13 , the intensification factor and the Nusselt number are compared between the 'Circular', 'Square' and circular 'Polished' mock-ups (defined in table 1)

all of them being made by AM. It shows that the mock-ups have the same thermal behavior, whether the cross-sections are circular or square; whether the surfaces are as-built or polished.

Moreover, the heat exchange increases with the flow rate, linearly in the double-logarithmic scale of the Nusselt number as a function of the Reynolds number. It is then possible to build the following correlation:

𝑁𝑢 = 0.20 × 𝑅𝑒 . × 𝑃𝑟 ( 16
)
Where the Pr number power is chosen equal to 1/3, to keep the same temperature dependency as the Colburn correlation.

The relative error of the correlation (eq. 16) when compared to our experimental data is between -16% and +10% and its coefficient of determination is R²=98%. Its validity domain is:

200 < 𝑅𝑒 < 6,600 (17)
2.6 < 𝑃𝑟 < 5.7

Hydraulics performances

The hydraulic performances of AM mock-ups are shown for the linear pressure drop (figure 14) and the Fanning number (figure 15) as a function of the flow velocity and the Reynolds number respectively, taking into account the geometrical factor uncertainty (right) or not (left).

When the geometrical factor uncertainty (right plots) is introduced, the comparisons can hardly be done, its effect (right plots) being very high compared to left plots differences. This first result confirms the importance of controlling the channel geometry and specially the hydraulic diameter during the manufacturing step [START_REF] Wildgoose | Impact of Additive Manufacturing on Internal Cooling Channels With Varying Diameters and Build Directions[END_REF], [START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF], [START_REF] Arie | Experimental characterization of an additively manufactured heat exchanger for dry cooling of power plants[END_REF][START_REF] Niknam | Additively manufactured heat exchangers: a review on opportunities and challenges[END_REF][START_REF] Da Silva | Thermal and hydrodynamic analysis of a compact heat exchanger produced by additive manufacturing[END_REF][START_REF] Dahmen | Characterization of channels made by laser powder bed fusion and binder jetting using X-ray CT and image analysis[END_REF][START_REF] Snyder | Understanding Laser Powder Bed Fusion Surface Roughness[END_REF]. Otherwise, as illustrated in figures 14 and 15 (right side), prediction of hydraulic performances seems compromised even for a 10% variation on the hydraulic diameter (±200 µm in the present study).

The mock-ups have been nevertheless manufactured with the same orientation on each baseplate.

Assuming that the geometrical factor shifts more or less all the curves in the same direction, it may be still possible to draw some general tendencies from the left plots. Looking at figure 14 and figure 15 from this viewpoint, the results may be discussed as follows.

Comparing 'Circular' and 'Circular bis' which differ by the printing parameters (table 1), a conclusion would be that the powder quality and the heating of the baseplate have a strong sensitivity on the friction loss (30%).

Comparing the 'Circular bis' and the 'Polished' mock-ups, for which powder and printing parameters are the same, provides information on the polishing effect. According to literature [START_REF] Da Silva | Thermal and hydrodynamic analysis of a compact heat exchanger produced by additive manufacturing[END_REF], [START_REF] Idel'cik | Mémento des pertes de charge[END_REF] polishing conducts to a lower Fanning number (around 15% for high Reynolds numbers when the roughness decreases from 20 to 2-3 µm, [START_REF] Da Silva | Thermal and hydrodynamic analysis of a compact heat exchanger produced by additive manufacturing[END_REF]) for the polished surface than for the as-built surface (by reducing the roughness). However, when comparing the pressure loss of the 'Circular bis' and the 'Polished' mock-ups, they are quite similar and the chemical polishing does not seem to have any influence on the friction factor. This may be due to the small roughness variation, in our case, between chemical polishing and raw channels with respective roughnesses of 6.7-7.7 µm and 4.1 µm.

Comparing the 'Circular' and the 'Square' mock-ups provide information on the crosssection geometry effect itself. At the same Reynolds number and at Re higher than 1,500, 'Circular' mock-up has a Fanning number 16% lower than the 'Square' mock-up. This result confirms the interest of changing the manufacturing process and using the freedom of geometry printing of the AM technique.

Global thermal-hydraulic comparisons

The thermal-hydraulic behaviors of the 'Circular', 'Square' or 'Polished' mock-ups are compared in figure 16. No significant difference is found between these three mock-ups. This means that both chemical polishing and cross-section have no influence on the intensification factor for a given pumping power.

Numerical study

This section describes the mock-ups Computational Fluid Dynamics (CFD) modeling setup and results. Ansys Fluent 2021R1 is used with Ansys Workbench to draw, mesh and calculate the studied case.

Computational domain geometry

In order to take into account the conduction effects, the computational domain geometry includes the solid part around the channel described in the 'description of the mock-ups' section.

The fluid computational domain geometry (figure 17) includes the channel itself to which upstream and downstream straight pipes have been added to numerically stabilize the flow. Looking at the CT-scans of the cross-section (figure 11), the channel volume can be seen as a sweep function of a half-circle (1.095 mm radius) and a partial circle (0.72 mm height). Based on this interpretation, the dissymmetry of the cross-section is taken into account in the fluid computational domain geometry (figure 18).

Numerical model set-up

Boundary conditions

On the external faces in contact with the cold fluid, a free stream temperature and a convective coefficient are set. The convective coefficient is calculated using the Colburn correlation

(23,670 W•K -1 •m -2 )
. At the inlet, a flow rate and a temperature are set. At the outlet, a pressure boundary condition is set with a prevention of backflow and a targeted flow rate specification.

Mesh

The meshing is done using Ansys Fluent Meshing. It is made of poly-hexcore cells for a total of 0.68Mcells in the fluid domain (figure 19). The meshing is refined at the walls in order to resolve the main flow phenomena and to avoid the wall function call, with non-dimensional distance values from the first fluid computational cell to the wall (y+) between 0.2 and 2•10 -3 .

A mesh sensitivity has been conducted for the simulation on the whole range of present Reynolds numbers (between 200 and 6000). The grids were made successively of 0. between the inlet and the outlet and the pressure drop varied both by less than +/-5%. We concluded that the resolution with 0.68 Mcells was sufficient for the present study.

Turbulence modeling

The Reynolds Averaged Navier-Stokes (RANS) ,k-𝜔 SST (Shear Stress Transport) model is used. The authors [START_REF] Mazzei | CFD Modelling Strategies for the Simulation of Roughness Effects on Friction and Heat Transfer in Additive Manufactured Components[END_REF] showed the accuracy of this model in the present context of both a roughness modeling approach as defined in the introduction and AM mock-up calculations. A turbulence model is needed because of the zigzag pattern which decreases the laminar/turbulent transition around Re=200 [START_REF] Anxionnaz-Minvielle | Influence of the meandering channel geometry on the thermo-hydraulic performances of an intensified heat exchanger/reactor[END_REF], [START_REF] Rosaguti | Laminar Flow and Heat Transfer in a Periodic Serpentine Channel[END_REF]. Thus, even if the simulations are done at low Reynolds numbers (Re between 200 and 6 000), the flow is assumed turbulent.

Roughness modeling

As shown by [START_REF] Stimpson | Roughness Effects on Flow and Heat Transfer for Additively Manufactured Channels[END_REF], the roughness value Ra depends on the surface orientation. However, the wall roughness is assumed to be constant through the domain in the present approach. It is a first source of uncertainty. Obviously a second source of uncertainty is the roughness measurement itself in similar conditions [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF]. When it comes to the present roughness modeling, a third source of uncertainty is the translation of roughness measurements into an ESGR concept. It is a subject of research since at least five decades ago [START_REF] Koch | Loss Sources and Magnitudes in Axial-Flow Compressors[END_REF][START_REF] Croner | Industrial use of equivalent sand grain height models for roughness modelling in turbomachinery[END_REF][START_REF] Adams | A Simple Algorithm to Relate Measured Surface Roughness to Equivalent Sand-grain Roughness[END_REF]. The option chosen in the present study, as a first approach, is to use the Stimpson correlation [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF] which was built in a similar AM context:

𝑘 𝐷 = 18 • 𝑅 𝐷 -0.05 (19)
However some discrepancies exist: the correlation was established in straight channels (wavy channel in the present study), printed with an angle of 45° (15° in the present study) for hydraulic diameters between 0.415mm and 1.246mm (2 mm in the present case).

In the present study, the Ra measurements of table 1 for the Circular section are used, which give a range from 6.7 to 7.7 μm. As previously mentioned, these measurements were done on the smoothest walls which were vertically manufactured and which represent only a small portion of the total walls.

Where 𝛽 * = 0.09 is a constant used in the SST k-𝜔 model. 𝑘 is the non-dimensional ESGR parameter (-) and 𝑦 the non-dimensional distance to the wall (-). It is defined as:

𝑦 = 𝜌 • 𝑢 • 𝑑 𝜇 (22) 𝑘 = 𝜌 • 𝑢 • 𝑘 𝜇 ( 23
)
Where 𝜌 is the density (kg•m -3 ), 𝑢 the friction velocity (m•s -1 ), 𝑑 the distance to the wall (m),

𝑘 the ESGR parameter (-) and 𝜇 the viscosity (Pa•s). Then the values of k and 𝜔 at the wall are determined:

𝑘 = 𝑘 • 𝑢 (24) 𝜔 = 𝜌 • 𝜔 • 𝑢 𝜇 (25)

Numerical versus experimental results

Validation of the present model

The CFD calculations are compared with experimental data with regard to the Nusselt number (figure 20) and the Fanning factor (figure 21) as a function of the Reynolds number. ± 20% deviation areas are indicated (red dotted lines) in order to better visualize the CFD versus experiment comparison. As far as CFD uncertainty quantifications are concerned, they generally rely on CFD results versus experiment comparisons and sensitivity calculations as an input.

Therefore present results can be considered as a first step towards this direction. However it is beyond the scope of the present paper to perform all the necessary calculations for an uncertainty quantification and to perform the relevant post-treatment of them in this view. It is an area of ongoing research [START_REF] Wenig | CFD Uncertainty Quantification using stochastic spectral methods-Exemplary application to a buoyancy-driven mixing process[END_REF]- [START_REF] Hassan | Uncertainty quantification (UQ) for CFD simulation of OECD-NEA cold leg mixing benchmark[END_REF] which is not mature enough to be used in the context of the present study. As a logical consequence, no uncertainty quantification are added to the numerical results in this paper. The experimental data have been collected in the 'circular' mock-up. The data have been reduced assuming the deformation as illustrated in figure 18. The equivalent hydraulic diameter is thus 1.94 mm instead of 2 mm.

The prediction of the Nusselt number is satisfactory (figure 20), the systematic overestimation giving an average relative error of 10%.

The agreement on the Fanning number is satisfactory (figure 21) with again an average relative error of 10%. The point at Re=897 seems out of the trend, it is due to a numerical convergence issue.

The two main modifications presented in this paper in comparison with a CFD approach of a conventionally manufactured heat exchanger are about the roughness modeling approach and the use of a different geometry than the CAD guideline. Neither approach includes tuning parameters.

The present roughness modeling relies on the ESGR as described in the previous section and then on the Ra measurement.

The change of geometry is a feature to underline because it means that the printing and the dimensional characterization of the considered object is required prior to CFD modeling. To the authors' knowledge, it is difficult to predict the geometry differences between a CAD geometry and the final AM printed component. This fact is a break with the usual modeling in which the design of an object is assumed to be known before manufacturing it.

In summary, the present results could not have been obtained from a blind calculation of a new design AM object because of the necessary knowledge for modeling of both Ra and the geometry change between CAD and reality. In view of this issue, a sensitivity study on both parameters is presented in the following sections.

Sensitivity to the ESGR

The ESGR used in the present modeling (results in the previous section) has been set up from the Stimpson correlation [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF] with ks=116.9 μm .

The sensitivity of ks on the flow is shown for the Nusselt number (figure 22) and the Fanning number (figure 23). Six different ks have been considered (ks=40; 80; 116.9; 160; 200; 240) at three different typical Reynolds numbers (Re=380; 730; 1,290). Each Reynolds value is not exactly the same as the typical target value because the mass flow rate is set in the simulation, not the Reynolds number. As the fluid physical properties depend on the temperature, the Reynolds number can slightly vary (max 1.7%).

These graphs show the relationship between the thermal-hydraulic prediction of the numerical model and the ESGR which is especially strong for high Reynolds number. The ESGR sensitivity increases with the Reynolds number for the Nusselt number (+36% on the studied range for Re=1300) and the fanning numbers (+26% on the studied range for Re=1300). It shows that neglecting the roughness or inputting an erroneous roughness may lead to significant discrepancies at high Reynolds numbers. At low Reynolds numbers, it may not be such a concern. Despite the various hypotheses that differ between the present study and the work of Stimpson et al. [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF], the calculated ESGR from eq. ( 19) proposed by the authors [START_REF] Stimpson | Scaling Roughness Effects on Pressure Loss and Heat Transfer of Additively Manufactured Channels[END_REF] seems well adapted to predict the thermal-hydraulic performances of the present meandering milli-channel made by AM without any needs in tuning parameter.

Sensitivity to the cross-section change

The sensitivity to the cross-section change on the thermal-hydraulic performances is shown for the Nusselt number (figure 24) and the Fanning number (figure 25) as a function of the Reynolds number. The objective is to compare the present model which is based on the actual geometry observed on the CT scans with the CAD cross-section, the latter being most likely used during blind CFD calculations e.g. for design optimization.

The cross-section affects both the Nusselt and the Fanning numbers. On the Nusselt number, the mean error done by simulations is of 9.9% considering the visualized cross-section and of 16.6

% considering the CAD cross-section. About the Fanning number, the choice of the cross-section is more complex and looking at the mean error is not sufficient (visualized cross-section: 8.6% or 6.5% not considering point at Re=897 and CAD cross-section: 4.1%). Using the visualized cross-section allows a better estimation especially at low Reynolds number (Re<400). Nevertheless, it creates higher calculation instability as for example at Re=897. However, the most realistic choice is to use the visualized cross-section in simulations, which leads to a better agreement with the experimental Fanning number curve.

Conclusion

In the present work, heat exchangers made by AM with mini-channels and meandering flows are studied both experimentally and numerically.

Experimentally, different mock-ups obtained either by HIP or by L-PBF AM, with different cross-sections (circular and square) and different surface roughnesses (raw and polished) have been characterized and compared.

On the one hand, the circular cross-section has a benefit of up to 16% on the Fanning number with no difference on the Nusselt number, comparing to square cross-section. On the other hand, a clear trend is found on the Nusselt number and the intensification factor versus the Reynolds number: moving from an HIP to AM leads to an increase by 80% on the Nusselt number, due to the roughness effects of the AM. Present Nusselt number versus Reynolds number experimental data are embedded in a proposed correlation within a precision range of +10/-16% with a coefficient of determination equal to 98%. These data therefore contribute to see AM as a great opportunity to heat exchange intensification.

Concerning the Fanning number, a gain up to 37% was measured between AM and HIP.

Amongst the various AM mock-ups, it was shown that the hydraulic diameter uncertainty based on available cross-sections CT-scans hindered a clear identification of differences.

A numerical approach to calculate the AM circular cross-section mock-ups tests is also presented. It is based on RANS k-𝜔 SST turbulence modeling and a roughness modeling of the ESGR type. Inputting in this model the measured roughness data on the actual mock-up is enough to get satisfactory thermal-hydraulic results, providing that the channel geometry is carefully computed using CT-scans of the actual mock-up and not using the original CAD. The numerical results are then in agreement with the experiment within errors comparing to experiment of 10% in average for the Nusselt number and within a maximum of 20% for the Fanning numbers.

Thus, there is still progress required to be able to perform blind CFD calculation and thermalhydraulic performances prediction of a new design AM object because of the necessary knowledge for modeling of both the roughness and the geometry change between CAD and reality. In the present case, the geometry modification introduced a shift of 6.7 % on the Nusselt number accuracy. Indeed the mean error is of 9.9% and 8.6% respectively on Nusselt and Fanning numbers.

The roughness is not an issue for Reynolds numbers below 400, but for Reynolds numbers around 1300, the Nusselt number could vary as much as 36% in our sensitivity calculation range.

Future works could be dedicated to the prediction of both roughness and dimension deviation from defined AM printing parameters to be able to perform a priori both modeling and design of thermal-hydraulic performances of the AM components. b Optical measurement using a focus variation microscope on vertically manufactured walls.

Figure 6

 6 Figure 6 and figure 7 display the thermal performances of the two mock-ups, intensification

  and figure 9), using the linear pressure drop ∆𝑃 𝐿 ⁄ (bar•m -1 ) and the Fanning number (-) respectively. To assess accurate Fanning numbers, the hydraulic tests (dedicated to pressure drop measurements) are performed in adiabatic conditions. This means that the data presented in figures 6 and 7 (non-adiabatic tests to assess thermal performances) do not come from the same experimental set of data than those presented in figures 8 and 9 (adiabatic tests). Since the temperatures and thus the physicochemical properties of the flow (specially the viscosity) differ between the adiabatic and non-adiabatic sets of experiments (around 20°C for hydraulic tests and 60°C for thermal tests), at a given flow velocity, the Reynolds numbers are different between figures 7 and 9.

  mock-up, separating the hot fluid and the cold fluid, m ESGR Equivalent Sand Grain Roughness f Fanning number, measuring the friction behavior of a flow h Convective heat transfer coefficient, W.m -2 .K -1 HIP Hot Isostatic Pressing k Turbulent energy of the k-ω model, m 2 •s -2 ks Equivalent Sand Grain Roughness parameter, m L Half length of the straight part of a channel pattern, rate of the k-ω model, s -

a

  Mechanical measurement using a direct roughness meter on vertically manufactured walls.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 1234678910 Figure 1 : Schematic diagram of the test rig

Figure 11 Figure 12 :Figure 13 :Figure 14 .

 11121314 Figure 11 Cross-sections measured by CT-scan, on 'Circular' mock-up, either in straight lengths

Figure 15 .

 15 Figure 15 . Fanning number (friction factor) compared between AM mock-ups, with (right) or

Figure 16 :

 16 Figure 16 : Intensification factor as a function of pumping power, with (right side) or without

Figure 17 :Figure 18 :Figure 19 :Figure 20 :

 17181920 Figure 17 : Computational domain geometry

Figure 21 :

 21 Figure 21 : CFD model calculation versus experimental data: Fanning number (friction factor).

Figure 22 :Figure 23 :Figure 24 :

 222324 Figure 22: Nusselt number as a function of the ESGR, for different Reynolds numbers

Figure 25 :

 25 Figure 25 : Cross-section sensitivity, Fanning number as a function of Reynolds number.

Fig 1 Figure 1 : 2 Figure 2 : 3 Figure 3 : 4 Figure 4 : 5 Figure 5 Fig 6 Figure 6 : 7 Figure 7 : 8 Figure 8 : 50 Fig 9 Figure 9 : 51 Fig 10 Figure 10 : 52 Fig 11 Figure 11 : 53 Fig 12 Figure 12 : 54 Fig 13 Figure 13 :Fig 14 Figure 14 :Fig 15 Figure 15 :Fig 16 Figure 16 :Fig 17 Figure 17 :Fig 18 Figure 18 :Fig 19 Figure 19 :Fig 20 Figure 20 :Fig 21 Figure 21 :Fig 22 Figure 22 :Fig 23 Figure 23 :Fig 24 Figure 24 :Fig 25 Figure 25 :

 11223344556677885099511010521111531212541313141415151616171718181919202021212222232324242525 Fig 1

  Stimpson et al. [6] measured Ra roughness of different direction-facing manufactured channel surfaces, and found in the case of Inconel (coupons L-1x-In and L-2x-In) a ratio of 1.75 between the vertical surface and the 45 deg downward-facing surface. Consequently it leads us to use in the present case an average roughness of 12 μm, yielding to ks = 120 µm.

	This ESGR is used in the High Roughness (Icing) model which is natively implemented in
	Ansys Fluent. The icing model [35] uses the following equations to compute the non-dimensional
	k and ω on a wall, named kw + and ωw + :
	𝑘 = max 0;	1 𝛽 *	𝑡𝑎𝑛ℎ	𝑘 30 𝑙𝑛(10) 𝑙𝑛	+ 1 -𝑡𝑎𝑛ℎ	𝑘 125	𝑡𝑎𝑛ℎ	𝑘 125	(20)
	𝜔 =	300 𝑘	tanh	15 4𝑘	+	191 𝑘	1 -𝑒

Table 1 :

 1 Mock-up's dimensional characteristics.

			Variable printing
	Cross-		parameters
	Name	Ra (μm) Rz (μm)	Remark
	section		Powder	Baseplate
			quality	temperature
		6.7-7.4 b ;
	Circular Circular
		7.7 a

Table 3 :

 3 measurements uncertainties

	Sensors	Uncertainty (including	Measurement range
		standard deviation and
		sensor calibration)
	Temperature	±0.06°C	12.7°C -77.9°C
	Flow rate (hot fluid)	±0.4%	1.2 kg/s -30.9 kg/s
	Flow rate (cold fluid)		1 104 L/h -1 483 L/h
	Differential pressure	±1.3%	3.0 mbar -299.2 mbar

Funding

This work was part of the Famergie project supported by the CARNOT Institute (French government grant).


Gian Paramo
email: gianparamo@ufl.edu

Arturo Bretas
email: arturo.suman-bretas@cnrs.fr

Sean Meyn
email: meyn@ece.ufl.edu

Microgrid Frequency Stability: A Proactive Scheme Based on Dynamic Predictions

Keywords: Phasor Measurement Units, Power System Protection, Underfrequency Load Shedding, Particle Filters

The dynamic nature of microgrids introduces challenges in the context of frequency stability. This work presents a framework where the future state of microgrid frequency is predicted and corrective actions are optimized. Predictions are generated through Bayesian filters leveraging synchronized data acquired via PMUs. Taking a proactive approach makes it possible to optimize corrective actions considering dynamic system conditions. Testing is conducted via Matlab simulations. The performance of the solution presented in this work is compared to traditional load-shedding schemes and to predictive solutions found in literature. The results indicate that the proposed framework outperforms both. Some of the advantages of this framework include a reduction in amount of load dropped during compensation, the use of adaptive parameters which eliminates the need to simulate contingency conditions, and dynamic uncertainty quantification provided the particle filter.

I. INTRODUCTION AND LITERATURE REVIEW

The flexibility and dynamic nature of microgrids pose challenges in terms of protection and control. Creating a protection scheme that can cover the many and rapidly changing operating conditions often observed in these systems is difficult [START_REF] Hooshyar | Microgrid Protection[END_REF]. Another important, yet unresolved, issue effecting the operation and reliability of microgrids is the need to find generation-load balance after an island is formed [START_REF] Farrokhabadi | Microgrid Stability Definitions, Analysis, and Examples[END_REF]. The process of attaining generation-load balance for a microgrid during islanded operation is fundamentally the same as finding generation-load balance in the macrogrid: Sources with available capacity are actuated and excess load is dropped [START_REF] Farrokhabadi | Microgrid Stability Definitions, Analysis, and Examples[END_REF]. The limitations of state-of-the-art frequency stability schemes are also the same at both, the micro and the macro level. Currently, the majority of automatic underfrequency load shedding (UFLS) solutions are built on decentralized architectures with little to no feedback, and operate at the feeder level [START_REF] Horowitz | Power System Relaying, 4th[END_REF]. Traditional UFLS schemes are reactive in the context that corrective actions are taken only after the system has entered a critical state. This leads to two commonly observed issues: delayed response and overshedding [START_REF] Sigrist | A review of the state of the art of UFLS schemes for isolated power systems[END_REF].

In order to overcome some of these challenges, the solution presented in [START_REF] Mynam | Islanding Detection and Adaptive Load Shedding[END_REF] utilizes an adaptive algorithm based on power measurements to drive a load-shedding process. In [START_REF] Javadi | Frequency Stability Constrained Microgrid Scheduling Considering Seamless Islanding[END_REF], a mixed integer linear programming (MILP) algorithm is used to optimize a cost function after a drop in frequency is detected via the rate of change of frequency (RoCoF). An adaptive, multi-stage load shedding controller is presented in [START_REF] Marzband | Adaptive load shedding scheme for frequency stability enhancement in microgrids[END_REF], where local voltage and frequency measurements are used to identify locations where compensation is required. In [START_REF] Kerdphol | Self-Adaptive Virtual Inertia Control-Based Fuzzy Logic to Improve Frequency Stability of Microgrid With High Renewable Penetration[END_REF], a fuzzy logic control algorithm is used to drive the virtual inertia in a microgrid in response to changes in power injections and system frequency oscillations. A cooperative control approach of neighboring distributed energy resources (DERs) is introduced in [START_REF] Trevizan | Cooperative Control of Energy Storage for Transient Stability Enhancement[END_REF]. In [START_REF] Gu | Cooperative Control to Enhance the Frequency Stability of Islanded Microgrids with DFIG-SMES[END_REF], neural networks are used to optimize the response of a cooperative control scheme for microgrid stability.

Predictive schemes in power systems have been suggested in various forms. A model predictive control approach is developed in [START_REF] Parisio | A Model Predictive Control Approach to Microgrid Operation Optimization[END_REF] to optimize microgrid operation. In [START_REF] Bretas | Real time instability prediction through adaptive time series coefficients[END_REF], measurements are collected via Phasor Measurement Units (PMUs). Data is processed as a time series, and short term predictions are made. The work presented in [START_REF] Larsson | Predictive frequency stability control based on wide-area phasor measurements[END_REF] predicts post disturbance steady state frequency. Active power injections of synchronous condensers are used by [START_REF] Sauhats | Enhancing Power System Frequency with a Novel Load Shedding Method Including Monitoring of Synchronous Condensers' Power Injections[END_REF] to detect loadgeneration imbalances proactively. Meanwhile [START_REF] Rudez | WAMS-Based Underfrequency Load Shedding With Short-Term Frequency Prediction[END_REF], develops a method where PMU measurements and polynomial curvefitting are used to make predictions.

In [START_REF] Pulendran | Energy Storage System Control for Prevention of Transient Under-Frequency Load Shedding[END_REF], frequency control is carried out by actuating fast acting DERs through an algorithm based on the Extended Kalman filter. In [START_REF] Farantatos | A Predictive Generator Out-of-Step Protection and Transient Stability Monitoring Scheme Enabled by a Distributed Dynamic State Estimator[END_REF], a predictive stability scheme was presented where the center of oscillations and the total energy of the system are used to detect unstable conditions shortly before the system enters a critical condition.

Solutions found in literature have the following limitations in common: Most cannot be supported by existing technology. For instance, these approaches normally rely on PMU measurements with high reporting rates (in some cases over 100 frames per second). This work, by contrast, is based on PMU reporting rates of 30 frames per second, as this rate would place a lighter burden on communication infrastructures [START_REF] Bose | Smart Transmission Grid Applications and Their Supporting Infrastructure[END_REF]- [START_REF] Naduvathuparambil | Communication delays in wide area measurement systems[END_REF].

Second, the techniques surveyed above fall into one of two categories: Data driven solutions and physics-based solutions. This work proposes a hybrid approach where data complements physics-based models. This enables the solution to produce results that take into account system dynamics and the physical limits of the system. The end product is a proactive framework for automatic frequency stability where model uncertainty is continuously updated [START_REF] Elfring | Particle Filters: A Hands-On Tutorial[END_REF].

The main contribution of this work lies in the application of the PF in the area of microgrid stability. The PF is leveraged to produce a hybrid solution that integrates a data-driven approach with physics-based equations.

The rest of the paper is structured as follows: Section II provides an overview of the technique, including an introduction to the particle filter. The complete scheme is summarized in Section III. Case studies and results can be found in Section IV. Concluding remarks are given in Section V.

II. DYNAMIC PREDICTIONS

A. Particle Filter

The particle filter (PF) is an approach to state estimation that combines Bayesian inference with the Monte Carlo method. The premise behind this type of estimation is that, it is easier to approximate a probability distribution, than it is to approximate an arbitrary nonlinear function or transformation [START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filters and estimators[END_REF]. Two advantages of the PF compared to other filtering techniques are that the PF is not limited to Gaussian distributions. The PF also offers inherit uncertainty quantification [START_REF] Paramo | Proactive Frequency Stability Scheme: A Distributed Framework Based on Particle Filters and Synchrophasors[END_REF]. The PF estimates an underlying probability distribution via samples called particles. Using a simple object tracking problem as an example, the process can be summarized as follows: First, a set of particles are placed around the predicted location of the object. Then measurements are used to assign weights to the particles. Particles closer to the measurements are given a higher weight. The object's location can be estimated by calculating the weighted average of the particles. The cluster of particles is then updated in two ways: First, particles with lower weights are redistributed around particles with higher weights. Then, the cluster itself is moved in accordance with the measurements and a system model. The process is then repeated at each time step. This is illustrated in fig. 1.

The problem of tracking the frequency of a power system is somewhat similar to tracking a particle moving cross a 2-D plane. For frequency tracking, the x-axis corresponds to time, while the y-axis corresponds to the frequency value at that particular moment. A sinusoidal function and the change in frequency, in this case ∆y, are part of the system model in this work. A complete derivation of the PF is omitted due to space constraints, but it can be found as part of the tutorial presented in [START_REF] Elfring | Particle Filters: A Hands-On Tutorial[END_REF]. For the interested reader, a comparison among several types of Bayesian filters can be found in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF].

B. Prediction

RoCoF is used as the primary means of disturbance detection. RoCoF offers a reasonable balance between adaptability and ease of implementation. After a disturbance is identified via the RoCoF, a prediction is made by feeding Artificial Data Points (ADPs) to the PF. ADPs are used to extend the horizon of the predictions produced by the PF beyond the typical k +1 time steps. This is an iterative process that creates a vector that is then used by the filter to make predictions:

ADP i = ADP i-1 + t s f ′ + t 2 s f ′′ (1)
ADP i-1 refers to the previous ADP. The average first derivative of frequency is represented by f ′ , while the average second derivative of frequency is represented by f ′′ . Finally, t s represents the time window of the derivatives [START_REF] Paramo | Proactive Frequency Stability Scheme: A Distributed Framework Based on Particle Filters and Synchrophasors[END_REF], [START_REF] Paramo | Proactive Frequency Stability Scheme via Bayesian Filters and Synchrophasors[END_REF].

C. Load Excess Calculation

Following a prediction, a load excess factor is calculated using a formula derived from the swing equation [START_REF] Horowitz | Power System Relaying, 4th[END_REF]:

L = R p H(1 - f 2 p f 2 1 ) p(f p -f 1) (2)
In which L represents the load excess factor, H is the inertia coefficient, p represents the power factor, and R p is the predicted rate of change of frequency.

R p = f p -f 1 t s (3)
During a prediction, R p is calculated as the difference between the current frequency measurement f 1 , and the predicted frequency measurement f p , divided by the time window t s . In this work, the state of frequency is predicted one second into the future.

III. SCHEME OVERVIEW

A. Distributed Operation

This works takes a distributed approach where areas (islands) are established based on PMU observability, as well as generation capacity. There is cooperation between neighboring areas, as PMUs near the edge of one area share information with adjacent areas.

PMUs are assumed to have a reporting rate of 30 measurements per second. Similar techniques in literate assume delays of roughly 0.5 seconds or less [START_REF] Larsson | Predictive frequency stability control based on wide-area phasor measurements[END_REF], [START_REF] Rudez | WAMS-Based Underfrequency Load Shedding With Short-Term Frequency Prediction[END_REF]. This delay accounts for communication latency and delays related to equipment operation, such as breaker trip times. This work makes the same assumption regarding the delay of the response. The local 979-8-3503-1360-4/24/$31.00 ©2024 IEEE inertia of each area is computed online as presented in [START_REF] Zeng | Online Estimation of Power System Inertia Constant Under Normal Operating Conditions[END_REF]. This value is then used to calculate the load excess factor via equation 2.

B. Optimized Response

A MILP algorithm is used to optimize the response after a disturbance. Three sources of compensation are considered in this work, these include: The actuation of DERs, shedding of non-critical loads (residential customers), and shedding of sensitive loads (industrial customers). Each one of them carries a different penalty for use, with DERs having the lowest penalty, and sensitive loads carrying the highest. The loading and capacity of the sources of compensation are updated in real time. This information is already captured by most utilities. After a suitable combination of compensation agents has been identified, trip commands are sent to the selected breakers. The process flow, including critical steps are presented in fig. 2.

IV. CASE STUDY

The case study presented in this section is based on a modified version of the IEEE-13 Bus test system used in [START_REF] Yu | Consensus Control of Distributed Energy Resources in a Multi-Bus Microgrid for Reactive Power Sharing and Voltage Control[END_REF]. The dynamics are modeled in Simulink, while the analytics and predictions are performed in Matlab. A One-Line diagram of the system is provided in fig. 3. System parameters, and perphase loading information can be found in tables 1 and 2 of [START_REF] Yu | Consensus Control of Distributed Energy Resources in a Multi-Bus Microgrid for Reactive Power Sharing and Voltage Control[END_REF]. Gaussian noise was added to the PMU measurements. It must be noted that the PF performs better when there is a moderate amount of noise in the measurements. This is because noise is related to how particles are distributed during the predict and update steps [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. To account for processing and communication delays, corrective actions are taken half a second after each prediction is made.

A. Case Study: Sub-islanding is Prevented

Initial conditions are as follows: All sources, including utility and DGs are online. The utility breaker opens at 1 second, creating an island with four DGs and a load excess of around 20%. Due to topology and physical distances, this larger island can be broken into two smaller islands as depicted in fig. 3. The scheme will attempt to keep smaller islands from being formed by keeping the DGs synchronized until frequency is stabilized. DG rotor angles will be tracked during this case study. Fig. 4 depicts the complete response of the system when the PF-based scheme is used to mitigate these disturbances. The prediction and mitigation process is as follows: An initial prediction is made after RoCoF thresholds are surpassed. This is illustrated in fig. 5. Half a second after the prediction is made (to account for communication and equipment delays), the first stage of compensation is carried out. Immediately after this stage of compensation is completed, a new prediction is made. This is shown in fig. 6. With frequency returning to normal ranges, the compensation algorithm stops. The complete response is shown in fig. 4, while difference in rotor angles is shown in fig. 10. The proposed scheme was able to stabilize the system while keeping all four DGs synchronized. In fig. 10,d-theta values are the rotor angle differences between DG1 and the other three machines.

B. Traditional UFLS Performance

A traditional UFLS scheme is used to mitigate the same disturbances. Fixed settings are used, and the complete response is shown in fig. 11. In this case the traditional UFLS scheme was able to stabilize frequency, however, synchronization between the DGs is lost as illustrated in fig. 12. In addition to maintaining synchronization between the islands, the solution presented in this paper produced corrective actions that are optimized taking into account dynamic system conditions. The traditional UFLS scheme on the other hand, simply drops load indiscriminately as frequency thresholds are surpassed [28]. In this work it was assumed that there was enough load available for shedding for the traditional UFLS scheme to stabilize the system.

C. Polynomial Curve Fitting

Finally, the framework presented in [START_REF] Rudez | WAMS-Based Underfrequency Load Shedding With Short-Term Frequency Prediction[END_REF] is used to mitigate the same disturbances. In [START_REF] Rudez | WAMS-Based Underfrequency Load Shedding With Short-Term Frequency Prediction[END_REF], predictions are produced via Polynomial Curve-Fitting (PCF). The limitations of the PCFbased solution become clear at the 2 second mark in fig. 13, and once again at the 4 second mark in fig. 14. The PCF-based solution is unable to leverage the dynamics of the system the way the PF-based solution can. As a result of this, in both cases, the PCF-based solution predicts that frequency will continue to decline and the scheme will 979-8-3503-1360-4/24/$31.00 ©2024 IEEE continue to compensate, leading to an unnecessary loss of customers. By contrast, as shown in fig. 7 and in fig. 9, the PF-based scheme is able to make predictions that decrease the amount of compensation utilized. In simulation based testing the method presented in this paper was able to stabilize the system while shedding 25% to 50% less load than the PCFbased technique. These results are consistent with the results obtained in [START_REF] Paramo | Proactive Frequency Stability Scheme: A Distributed Framework Based on Particle Filters and Synchrophasors[END_REF].

V. CONCLUSION

A proactive frequency stability framework for microgrids was presented in this work. After a disturbance was detected the future state of frequency was predicted via particle filters. Being able to predict a critical state a few seconds in advance made it possible to optimize corrective actions. During testing this solution outperformed UFLS schemes, as well as frameworks found in literature by reducing the amount of load shed in response to disturbances. Future work will focus on validation on larger and more complex systems.

Fig. 1 :

 1 Fig. 1: Particle Based Sampling.

Fig. 2 :

 2 Fig. 2: Process flow.

Fig. 3 :

 3 Fig. 3: One Line Diagram of the Test System.

Fig. 4 :

 4 Fig. 4: Frequency Response with PF Mitigation.

Fig. 5 :

 5 Fig. 5: Initial PF Prediction.

Fig. 6 :

 6 Fig. 6: Second PF Prediction.Since this new prediction shows that frequency is still on the decline, a new stage of compensation and prediction is carried out at around 2 seconds. This is illustrated in fig.7.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Third PF Prediction.This time however, the prediction shows the system's frequency climbing back to normal ranges. A subsequent stage of compensation is put on hold until the a new decline in frequency is detected at roughly 3.5 seconds, as shown in fig.8.

Fig. 9 :

 9 Fig. 9: Final PF Prediction.

Fig. 10 :

 10 Fig. 10: Rotor Angles with PF-based Mitigation.

Fig. 11 :

 11 Fig. 11: Frequency Response with Traditional UFLS Scheme.

Fig. 12 :

 12 Fig. 12: Rotor Angles with Traditional UFLS Mitigation.

Fig. 13 :

 13 Fig. 13: PCF-based Prediction.

Fig. 14 :

 14 Fig. 14: PCF-based Prediction.

This work was partially funded by the U.S. Department of Energy under Contract DE-AC05-76RL01830.


Anthony Chambon

Zoé Anxionnaz-Minvielle

Grégory Cwicklinski

Nathalie Guintrand

Aline Buffet

Bernard Vinet

SHELL-AND-TUBE HEAT EXCHANGER GEOMETRY MODIFICATION: AN EFFICIENT WAY TO MITIGATE FOULING

Crude oil fouling of a shell-and-tube heat exchanger sized according to TEMA standard is compared to a No-Foul design under industrial operating conditions. For similar operating conditions, TEMA and No-Foul heat exchangers have the same behavior regarding fouling. Since the No-Foul one has less tubes by design for the same heat duty, shear stress is increased. Consequently, the No-Foul heat exchanger is less prone to fouling at the same throughput. Impact of tube bundle geometry is then investigated. Helicallyfinned tubes are compared to plain tubes in the No-Foul heat exchanger. Under similar operating conditions, fouling rates measured are up to an order of magnitude lower than plain tubes (respectively 10 -11 and 10 -10 m 2 K/J). However, pressure drop across the tube-side in both No-Foul plain and finned setup are increased in comparison to the TEMA heat-exchanger.

INTRODUCTION

Fouling of heat exchangers is the build-up of fouling layers on the heat transfer surface. Crude oil fouling in refineries causes several operating, financial and safety issues. Before entering the atmospheric column, crude oil is usually heated in the preheat train by recovering heat from hot streams. But thermal performance is reduced by fouling due to an increase of the thermal resistance leading to additional fuel consumption and CO2 emissions. Moreover, the gradually decrease in tube cross-sectional area due to the fouling layer growth requires more pumping power to maintain the throughput. Ultimately fouling can cause a fluid blockage. Before reaching this stalemate, cleaning of heat exchangers is required. It could be performed using chemical or by dismantling and water-jetting the equipment. All these consequences are the root of numerous additional expenditures. Costs linked to fouling are estimated to USD3.6 billion a year [START_REF] Coletti | Crude oil fouling, deposit characterization measurements and modelling[END_REF].

Unfortunately, fouling phenomenon remains still not well understood. Therefore, several mitigations techniques have been explored to manage fouling in preheat train such as optimization of operating conditions, preheat train monitoring to optimize cleaning cycles [START_REF] Müller-Steinhagen | Heat exchanger fouling: mitigation and cleaning strategies[END_REF][START_REF] Ishiyama | Scheduling cleaning in crude oil preheat train subject to fouling: incorporating desalter control[END_REF], use of chemicals [START_REF] Baxter | Dispersant material for mitigating crude oil fouling of process equipment and method for using same[END_REF][START_REF] Brant | Polyalkyl succinic anhydride derivatives as additives for fouling mitigation in petroleum refinery processes[END_REF], tube inserts [START_REF] Aquino | Towards an understanding of how tube inserts mitigate fouling in heat exchangers[END_REF][START_REF] Krueger | Heat exchanger performance enhancement through the use of tube inserts in refineries and chemical plants -successful applications: Spirelf®, Turbotal®, and Fixotal® systems[END_REF]. This study focuses particularly on the technique of enhanced heat exchanger geometries.

Shell-and-tube heat exchangers are commonly set up in hydrocarbon services since it is perceived in industry as a safer option [START_REF] Coletti | Crude oil fouling, deposit characterization measurements and modelling[END_REF]. Nevertheless, fouling of other heat exchanger technologies have also been investigated so far. Compact and spiral heat exchangers successfully mitigate fouling [START_REF] Tamakloe | Design of Compabloc exchangers to mitigate refinery fouling[END_REF][START_REF] Wilhelmsson | Consider spiral heat exchangers for fouling applications[END_REF] by providing higher shear stress at a given flowrate. Shell-and-tube heat exchanger geometry modification (helical baffle, EM-Baffle or twisted tubes) can also help to mitigate fouling but improvement is limited to the shell-side where dead zones are suppressed [START_REF] Master | Fouling mitigation using Helixchanger ® heat exchangers[END_REF][START_REF] Brignone | EMBafle® in refinery service: on-field study and data validation through SmartPM®[END_REF].

However, tube-side geometry modification has received little attention up to date. Yang et al. [START_REF] Yang | Mitigation of crude oil fouling by turbulence enhancement in a batch stirred cell[END_REF] and Crittenden et al. [START_REF] Crittenden | Crystallization fouling with enhanced heat transfer surfaces[END_REF] studied the fouling of an enhanced surface (by the mean of an incorporated wire) and compared it to a smooth one. For the same operating conditions, fouling resistance measured were lowered with the enhanced surface. Provost et al. [START_REF] Provost | Field tests on heat exchangers equipped with dual enhanced tubes in a quench water application[END_REF] and Sippel et al. [START_REF] Sippel | Update on field tests on heat exchangers equipped with dual enhanced tubes in a quench water application[END_REF] published a feedback on fouling of internally helically finned tube with quench water of an industrial naphtha cracker. Enhanced tubes did not undergo any significant fouling after 24 months of unit operation. To the authors' knowledge, no fouling study of these tubes has been carried out with crude oil.

That's why we will study the influence of the sizing method and tube technology on fouling. Hence, fouling propensity of a shell-and-tube heat exchanger sized according to TEMA standards [START_REF] Tema | Standards of the Tubular Exchangers Manufacturers Association[END_REF] is compared to a heat exchanger sized with a No-Foul method as described by Bennett and Nesta [START_REF] Bennett | Reduce fouling in shell-and-tube heat exchangers[END_REF]. Impact of tube bundle technology is secondly addressed by testing successively plain and helically-finned tubes in the No-Foul heat exchanger.

The helically-finned tubes geometry features are illustrated in Fig. 1 [START_REF] Ji | Prediction of fully developed turbulent heat transfer of internal ribbed tubes -an extension of Gnielinski equation[END_REF]. Compared to the plain tubes, internal area is increased by 30% (i.e. Ai,fin / Ai,plain = 1.30).

EXPERIMENTAL SETUP

Fouling rig

The heat exchanger to be tested (the test section) is settled in a fouling rig described successively by Ratel et al. [START_REF] Ratel | Investigation of fouling rates in a heat exchanger using an innovative fouling rig[END_REF] and Chambon et al. [START_REF] Chambon | Fouling in the crude oil distillation preheat train: comparison of experimental data with model results[END_REF] during past issues of this meeting. So, the features of the rig are not recalled here and the reader is strongly invited to refer to these former papers to learn about it.

As in the industrial process, atmospheric tower bottom (ATB) flows through the shell side to heat up the crude. Oil used for experiments was sourced from the Black Sea area. Physical properties and chemical composition of both fluids are respectively reported in Tables 1 and2.

Crude oil and ATB used for No-Foul tests come from a different batch than that used for TEMA tests, nevertheless their composition shows a very strong similarity. However, the ATB used for the No-Foul runs would appear to contain less sulfur and asphaltenes which are known as fouling precursors. As a consequence, the ATB should have, a priori, a less fouling propensity than the one used for TEMA tests.

Heat exchanger designs

The sizing method originally proposed by the TEMA is known to promote fouling. Fouling resistances listed by the TEMA are fixed values whereas the fouling is a dynamic phenomenon which depends on the operating conditions. These inappropriate values lead to oversize heat exchanger which, finally, amplifies their fouling. Indeed, the addition of extra heat transfer surface tends to reduce flow velocities and the convective heat transfer coefficients and finally increase the wall temperatures. That's why Bennet and Nesta (2004) have developed a sizing method that takes better account of fouling.

The No-Foul test section is sized to fulfill their recommendations:

 tube-side wall temperature below 300°C and flow velocity higher than 2 m/s at rated throughput;

 on the shell side, leakage currents are minimized and cross-flow (B-stream) promoted;

 over-sizing is reduced to a 15% extra heat transfer area and no fouling resistance is considered in the design.

A comparative datasheet in table 3 gives the geometrical characteristics of the exchangers tested. Each tubular heat exchanger tested is made with carbon steel and has a similar heat flux level and the same tube gauge than the heat exchanger of a refinery pre-heat train.

Consequently, at constant flowrate, the No-Foul heat exchanger requires less heat transfer area than the TEMA in order to transfer the same duty. This results in fewer tubes. Hence, the flow velocity and the wall shear stress are increased, promoting suppression of fouling deposit. However, the No-Foul exchanger will generate more pressure drop than the TEMA due to the higher flow velocities.

Operating conditions range

Crude oil velocities and film temperature investigated respectively have been selected to be representative of industrial operating conditions.

Tests are performed in order to investigate effect of film temperature and velocity separately. After each fouling run, both sides of the test heat exchanger are flushed with gasoil for at least two days in order to remove as fouling material as possible.

Data reduction

Fouling resistances and fouling rates are derived from measurements of the overall heat transfer coefficient U(t). The latter is figured out with Eq. ( 1) by continuous monitoring of inlet and outlet temperatures. More details are given in the previous study [START_REF] Chambon | Fouling in the crude oil distillation preheat train: comparison of experimental data with model results[END_REF].

𝑄 = 𝐹 𝑈(𝑡) 𝐴 𝐿𝑀𝑇𝐷 (1)
In Eq (1), Q is the mean duty derived from hot and cold-side heat balances. LMTD is the logarithmic mean temperature difference and F the correction factor. In order to make reliable fouling propensity comparisons fouling rate and resistances must be related to the same heat transfer surface. Since plain and helically-finned tubes have the same outer diameter (1 inch), the external heat transfer area (Ao) is chosen.

Hence, all the fouling resistance is affected outside the tubes albeit it occurs on both sides of the heat exchanger.

Tube wall shear stress, τ, is figured out with Eq. ( 2), where ρ and u are respectively crude oil density and velocity.

𝜏 = 𝑓𝜌𝑢 (2)
The Fanning friction factor f is respectively calculated with Blasius equation (Eq. 3) for plain tubes and with the Zdaniuk et al. [START_REF] Zdaniuk | Experimental determination of heat transfer and friction in helically-finned tubes[END_REF] correlation (Eq. ( 4)) for helically-finned tubes. In these equations, Re refers to the crude Reynolds number and N, e, α respectively to the fin height, the number of fin starts and the helix angle (see Fig. 1).

𝑓 = 0.0791𝑅𝑒 . (3)
𝑓 = 0.128𝑅𝑒 . 𝑁 . (𝑒 𝑑 ⁄ ) . 𝛼 . (4)
The Zdaniuk et al. [START_REF] Zdaniuk | Experimental determination of heat transfer and friction in helically-finned tubes[END_REF] correlation has been chosen since it gives the better agreement with experimental pressure drop measurements made during a preliminary study.

According to Ebert and Panchal [START_REF] Ebert | Analysis of Exxon crude-oil-slip stream coking data[END_REF], the film temperature (Tf) is a weighted average between wall (Tw)

and bulk (Tb) temperature (Eq. 5).

𝑇 = 𝑇 + 0.55(𝑇 -𝑇 ) (5)
Bulk temperature is assumed to be the average between crude inlet and outlet temperature during the whole run duration. Wall temperature is derived from Eq. [START_REF] Aquino | Towards an understanding of how tube inserts mitigate fouling in heat exchangers[END_REF].

𝑇 = 𝑇 + (6)
Where hi is the internal convective heat transfer coefficient and Ai, the effective internal heat transfer area. Convective heat transfer coefficient is respectively figured out with Dittus-Boelter (Eq. 7) and Zdaniuk et al. [START_REF] Zdaniuk | Experimental determination of heat transfer and friction in helically-finned tubes[END_REF] correlation (Eq. 8) for plain and helically-finned tubes.

= 0.023𝑅𝑒 . 𝑃𝑟 .

= 0.029𝑅𝑒 . 𝑁 . (𝑒 𝑑 ⁄ ) . 𝑃𝑟 .

In these equations, k is the thermal conductivity of crude oil (see Table 1) and Pr, the crude Prandtl number. The fouling rate of tests NF-P2 and T-P7 are comparable. Although on NF-P2 the film temperature is higher, the crude velocity (wall shear stress) is lower and the ATB cross-flow velocity is 20% lower. The gap could be explained by a little bit lower fouling propensity of the ATB (see Table 2).

RESULTS

The absence of fouling for the NF-P3 and T-P2 runs proves that they are comparable. On test T-P2 a fouling resistance close to uncertainty was detected while no fouling was detected on test NF-P3. These observations seem to be consistent with the lower crude oil temperature and with the lower fouling propensity of the ATB.

Although test NF-P4 is carried out at a higher temperature than test T-P1, no fouling is noticed. The threshold temperature was probably not reached. Referring to Fig. 2, the real threshold curve would be positioned above NF-P4.

For a similar wall shear stress (i.e. oil velocity) on tests NF-P5 and T-P5, a lower fouling rate was observed for the NF-P5 test. This is in accordance with the lower film temperature and highlights the strong impact of temperature on fouling growth. It is confirmed by comparing NF-P5 and T-P6 (carried out at lower temperature and similar shear stress). Therefore, the real threshold curve is assumed to be located between these two points which is in accordance with the Panchal et al. [START_REF] Panchal | Threshold conditions for crude oil fouling, in Understanding Heat Exchanger Fouling and its Mitigation[END_REF] forecast (see Fig. 2).

At low flow velocities and film temperatures, the equilibrium between deposition and suppression rates is slower to be reached (slow dynamics). Conversely, at higher velocities and film temperatures, the equilibrium of the two mechanisms is more rapidly set up (fast dynamics). Hence, the Panchal et al. [START_REF] Panchal | Threshold conditions for crude oil fouling, in Understanding Heat Exchanger Fouling and its Mitigation[END_REF] correlation delineates the fouling and non-fouling areas poorly when the balance between the deposition and suppression rate is more difficult to set up. Forecasts are more accurate for high flow velocities when the dynamics is more pronounced.

To sum it up, for similar operating conditions (velocity and film temperature), the TEMA and No-Foul heat exchangers have the same behavior regarding fouling. Fouling rate and resistance measured are all included in measurement uncertainties.

Comparison at same throughput. Similarly to Fig. 2, results are gathered in Fig. 3 where horizontal axis is now the crude oil flowrate. This graph allows to compare the fouling propensity of a shell-and-tube heat exchanger according to the method used for its sizing.

The two threshold curves are plotted with the same correlation (Panchal et al. [START_REF] Panchal | Threshold conditions for crude oil fouling, in Understanding Heat Exchanger Fouling and its Mitigation[END_REF]) and the same optimized parameters. They are not superimposed because the threshold condition depends on the wall shear stress but not directly to the crude oil flow rate. Indeed, for the same flowrate, wall shear stress is greater in the No-Foul heat exchanger because crude oil velocity is higher due to fewer tubes. Consequently, the No-Foul heat exchanger threshold curve is translated horizontally to the left-hand side of the graph.

By replacing the TEMA heat exchanger with a No-Foul one, fouling operating loci located between the two thresholds curves are moved into the non-fouling area. Thus, fouling is reduced because operating conditions are now in the non-fouling zone of the dashed curve.

Cost considerations. No-Foul heat exchangers prove to be more economical to produce because they are sized as tightly as possible. Since both exchangers are tubular heat exchangers, replacing a TEMA with a No-Foul avoid heavy modification of the unit because the bulk is similar and nozzle locations can be kept in place. Finally, operating costs are reduced since the No-Foul is less prone to fouling. The main drawback of the latter design is pressure drop due to the higher flow velocities at same throughput.

Helically finned vs. plain tube

In this section, the plain tube bundle of the No-Foul heat exchanger was replaced by a bundle of helically finned tubes. The two bundles have the same features: same tube count, same baffle spacing, same outer diameter, etc. The tube internal geometry is the only modified parameter. The goal is to compare the fouling propensity of these tubes with respect to the plain tubes at similar operating conditions.

Results of fouling tests are summarized as a function of the wall shear stress in Fig. 4. Fouling runs are named with the 'NF-F' label. Similarly, fouling rates are expressed into brackets relatively to the T-P4 value. For readability reasons, only T-P4 and T-P8 are plotted for the TEMA-series.

Fouling was detected during NF-P2 run carried out with plain tubes while no fouling is measured at the end of the test NF-F1. Both tests are performed at a similar wall shear stress. Nevertheless, shear stress figured out with the Zdaniuk et al. [START_REF] Zdaniuk | Experimental determination of heat transfer and friction in helically-finned tubes[END_REF] correlation (Eq. 5) is derived from experimental pressure drop measurements which include losses caused both by friction and drag. Thus, the shear stress assessed is an overall wall shear stress. Overall wall shear stress cannot explain the absence of fouling at the end of NF-F1 test. This suggests that local (rather than overall) wall shear stress acts in helically-finned tubes to prevent fouling.

The fouling rate and resistances measured on test NF-F2 are an order of magnitude lower than for the plain tubes (T-P4). Temperature monitoring of tests NF-F2 revealed that the tube-side was fouled first.

However, it cannot be known whether the shell-side was subsequently involved in fouling. The ATB processed on test T-P4 is assumed to be more prone to foul and could explain why fouling is more severe in this test but it cannot account for the large difference in fouling rates. The previous hypothesis of higher local shear stress in helically-finned is believed to explain fouling mitigation.

No fouling resistance was detected on the NF-F3 test, while a fouling resistance near of the uncertainty Helically-finned tubes would appear to be less prone to fouling than plain tubes, at least for the range of operating conditions investigated. Fouling of helically-finned tubes has only been noticed for very low velocities and high film temperatures. For similar operating conditions, the fouling rate measured are is one order of magnitude lower than for plain tubes. However, the relatively short duration of the tests (300 to 500 h) and the limited number of experimental runs do not allow to conclude about fouling over the long periods of time and on other operating ranges.

The comparison of the tests carried out under similar operating conditions suggests a different rate of foulant growth which could explain the lower fouling propensity of helically-finned tubes. Fins could generate boundary layer detachments and impart swirling flow motion. Thus, at same film temperature and shear stress, the enhanced local turbulence would reduce the fouling rate.

However, the mitigation effect of finned tubes is counterbalanced by a significant increase in pressure drop compared to the plain tubes. At similar flowrate, the higher shear stress causes more pressure drop.

CONCLUSIONS

1. For similar operating conditions (velocity and film temperature), shell-and-tube heat exchangers sized according TEMA standard and with a No-Foul design method have the same behavior regarding fouling.

2. At the same throughput, the No-Foul sized heat exchanger has proven to be effective to mitigate fouling.

3. Helically-finned tubes have a lower fouling propensity than plain tubes. Fouling rate measured are an order of magnitude lower than for the plain tubes (10 -11 and 10 -10 m 2 K/J respectively). However, the relatively short duration of the tests (300 to 500 h) and the limited number of experimental runs do not allow to conclude about fouling over the long periods of time and on other operating ranges.

4. The lower fouling propensity of helically-finned tubes could perhaps be explained by a higher local wall shear stress due to local turbulence enhancement.

5. At the same throughput, pressure drop is higher in the No-Foul sized heat exchanger than in the TEMA sized one. Similarly, our finned tubes bundle causes significantly more pressure drop than in the plain tube bundle. Fouling for sizing m 2 K/W 1.9 10 -3 0 Fig. 1 Geometry description of a helically-finned tube [START_REF] Ji | Prediction of fully developed turbulent heat transfer of internal ribbed tubes -an extension of Gnielinski equation[END_REF]. Fig. 4 Fouling experiment results for plain and helically finned tubes.

Fig. 1 Geometry description of a helically-finned tube [START_REF] Ji | Prediction of fully developed turbulent heat transfer of internal ribbed tubes -an extension of Gnielinski equation[END_REF].

  Fig. 2. Tests performed with the TEMA heat exchanger are labelled with 'T-P' whereas 'NF-P' refers to

  was measured on the T-P8 test. It is unlikely fouling on the T-P8 test is due to the ATB because under similar shell-side conditions, no fouling was measured in the T-P2 and T-P6 tests. Thus, the fouling resistance on the T-P8 test could only be due to crude oil. Since the NF-F3 and T-P8 tests are carried out at the same film temperatures and shear stress, the higher local shear stress seems, once again, to be the most likely hypothesis to explain the fouling mitigation.

6 .

 6 Tests carried out with a non-fouling fluid (alternately on the shell-side and on the tube-side) and a hydrocarbon fluid (crude or ATB) would allow to dissociate the contribution of each fluid in the fouling of the heat exchanger. Moreover, a better knowledge of mechanisms involved in crude oil fouling would provide additional evidences for understanding the lower fouling propensity of the helically-finned tubes.

Fig. 2 Fig. 3

 23 Fig. 2 Fouling experiment results for the TEMA and No-Foul heat exchangers as a function of crude oil

Fig. 2 Fig. 3 Fig. 4

 234 Fig. 2 Fouling experiment results for the TEMA and No-Foul heat exchangers as a function of crude oil

Table 1 .

 1 Physical properties of the oil and the ATB.

	Physical properties at 300°C

Table 2 .

 2 Chemical composition of the oil and the ATB.

	Fluids used for testing:	TEMA			No-Foul
		Oil	ATB	Oil	ATB
	Asphaltenes (%w)	1.1	2.5	1.1	0.46
	CCR (Conradson Residue)	3.92	5.96	4.45	5.19
	Sulphur (%w)	1.22	2.25	1.37	0.56
	Coking factor	0.81	2.05	0.77	1.25
	Cracking factor (%w)	1.92	2.20	1.99	2.37
	Σ C aromatics (%w/w)	17.5	24.5	16.2	17.5
	Σ C saturates (%w/w)	82.5	75.6	83.8	82.5

Table 3 .

 3 Geometrical characteristics of the shell-and-tube heat exchangers compared.

	Heat exchanger design		TEMA	No-Foul
	Tube bundle		Plain	Plain	Finned
	Tube length	m	2.0	2.0
	Tube outer diameter (do)	mm	25.4	25.4
	Tube inner diameter (di)	mm	19.9	19.9	21.4
	Tube count	-	32	28
	Tube pitch	mm	31.75	31.75
	Tube arrangement	-	square	square
	Tube pass	-	4	4
	Shell	-	1	1
	Shell internal diameter	mm	260	248
	Number of baffle	-	17	13
	Baffle spacing	mm	100	115
	Baffle window	%	20	24
	Tube-to-baffle clearance	mm	0.9	0.4
	Shell-to-baffle clearance	mm	2.0	0.9
	Heat transfer area (Ao)	m 2	5.1	4.3
	Heat duty	W	5 10 4	5 10 4
	Heat flux	W/m 2	10 4	1.2 10 4