Expanded Microchannel Heat Exchanger: Finite Difference Modeling

: A ﬁnite difference model of a heat exchanger (HX) considered maldistribution, axial conduction, heat leak, and the edge effect, all of which are needed to model a high effectiveness HX. An HX prototype was developed, and channel height data were obtained using a computerized tomography (CT) scan from previous work along with experimental results. This study used the core geometry data to model results with the ﬁnite difference model, and compared the modeled and experimental results to help improve the expanded microchannel HX (EMHX) prototype design. The root mean square (RMS) error was 3.8%. Manifold geometries were not put into the model because the data were not available, so impacts of the manifold were investigated by varying the temperature conditions at the inlet and exit of the core. Previous studies have not considered the inﬂuence of heat transfer in the manifold on the HX effectiveness when maldistribution is present. With no ﬂow maldistribution, manifold heat transfer increases overall effectiveness roughly as would be expected by the greater heat transfer area in the manifolds. Manifold heat transfer coupled with ﬂow maldistribution for the prototype, however, causes a decrease in the effectiveness at high ﬂow rate, and an increase in effectiveness at low ﬂow rate.


Introduction
Heat exchangers (HXs) are ubiquitous in modern life [1] and have a market size predicted to reach over $20 billion U.S. dollars by 2024 [2]. There has been substantial interest in improving heat exchanger performance by utilizing microchannel-based HX [3][4][5][6]. Even with improved performance, minimizing economic costs are a core driver of specific applications. For example, an emerging application, solar water pasteurization, demands low-cost HX and has the potential to save the lives of thousands of children a day [7][8][9]. Another important application can be found in absorption cycles, which can utilize waste heat from diesel engines, such as those found on ships [10].
One of the ways gaining prominence for reducing the costs of products is applying the open-source paradigm of technical development to hardware [11][12][13][14]. An open-source polymer laser welder has been developed [15], heat transfer in polymer HX using it was studied [16] and recent work [17] has also provided a low-cost method of making a polymerbased expanded microchannel HX (EMHX). EMHX modeling is more complicated than conventional microchannel HX [18]. A study performed a computed tomography (CT) scan to examine a prototype polymer-based EMHX such as this, and it indicates significant flow maldistribution due to channel height variations [19]. There have been significant improvements in the modeling of micro-channel HX, which can be used to drive the development of EMHX [4,[20][21][22][23][24][25][26][27][28][29][30][31]. Recently, a paper [32] developed a finite difference model that can be applied to this EMHX that was scanned, allowing for a comparison between modeled and experimental heat transfer characteristics when maldistribution is present.
The performance of a HX is characterized by its effectiveness, η, defined as the ratio of the actual heat transfer rate to the maximum heat transfer rate. (1) Here we use effectiveness to mean the approach to the ideal temperature, which can be applied to the HX as a whole, or individual channels. The actual heat transfer rate is: where the heat capacity rates, C H and C C , are defined using mass flow rate and specific heat.
m H C p,H , C C = . m C C p,C .
The maximum heat transfer rate is: where C min is the minimum heat capacity rate between the hot and cold side. The number of transfer units or NTU is defined as the ratio of the heat transfer rate between the fluids to the minimum heat capacity rate. NTU represents the heat transfer ability of a heat exchanger normalized by heat transfer area. A high NTU value indicates a low flow rate. NTU = U A C min .
Here A is the heat transfer area and U is the overall heat transfer coefficient. If channels in an HX vary in size, then fluid velocity will be greater in some channels. The channels with greater fluid velocity will have a smaller NTU, which means the fluid in these channels will not cool down as much. The mechanisms producing channel irregularity include nonuniform fouling [33], two phase flow [34], thermally induced maldistribution caused by unequal thermal expansion coefficients [35], and non-uniform spacing of fins [36].
To understand the way that fluid travels through the channels of the HX, a few definitions are required. Consider a cross section of the core of the HX in Figure 1.
When the flow is uniform, the flow rate is the same for each channel (hot and cold alike because it is balanced flow). Now, consider the cold and hot channels of a cross section as separate groups. Positive correlation is defined as a group of cold and hot channels having the same flow maldistribution, e.g., the right-side channels have lower flow rate for both cold and hot. Negative or anti-correlation is defined by hot channels having the opposite flow rate trend as the cold channels, e.g., the cold channels on the right side having lower flow rate than the cold channels on the left, but the hot channels on the left side having lower flow rate than the hot channels on the right. Uncorrelated flow occurs when the flow rates of the hot and cold channels show neither positive nor negative correlation, which could mean that one temperature of channels has constant flow rate and the other temperature of channels could have any pattern. Microchannels (<1 mm hydraulic diameter) that are smooth and contain a liquid with negligible viscous heat generation have the same physics as traditional laminar correlations [37][38][39][40][41]. Therefore, traditional laminar correlations are used for this paper.
Axial conduction is the conduction of heat in the direction of fluid flow, by both the fluid and the HX material [42,43]. This can be a significant source of ineffectiveness for high effectiveness HXs. The ineffectiveness due to axial conduction can be added as post-processing for a large range of configurations [44].
The heat leak is the gain or loss of heat from the environment. As a first approximation, the heat leak ineffectiveness can be added as a linear term to the other ineffectivenesses [45].
Edge effects are considered in that the channels near the edge have fewer neighbors and thus do not transfer as much heat, giving those channels a lower NTU. This edge effect actually penetrates to channels not on the exterior as well and can be important even for 80 layers [46].
Large temperature variations in the HX can significantly change the properties of the fluid and the wall [47]. However, the prototype HX for this paper has small temperature variations, so these effects are not considered.
A number of papers discussed flow maldistribution due to the manifold [48][49][50][51]. Different paths of the fluid in the manifold (see Figure 1) have different flow resistances, so channels have varying flow rates in the HX core. The problem of manifold heat transfer not being the same for all fluid paths has not been analyzed [52]. The present paper addresses this analysis by varying the boundary conditions at the HX core based on different manifold heat transfer characteristics to examine the effect on each fluid path in the core. With improvement to the model using the conclusions of the analysis, better fit between the modeled and experimental results can be achieved and design of microchannel heat exchangers can be improved by taking into account manifold heat transfer with maldistribution.

The Advanced Model
The advanced model is a finite difference model that calculates the flow in each channel of the HX based on the height of each channel as a function of its length. The advanced model takes into account edge effects, heat transfer between all neighboring channels, axial conduction, and heat leak. Edge effects and heat transfer between neighboring channels are considered with a transverse heat flux term, see Figure 2 and Equation (6): where A HT is the total heat transfer area shared with one neighbor, U is the overall heat transfer coefficient, and the temperature difference between the example cold channel and a transverse neighbor k is T tran,k − T c,i,j . While two of the transverse channel neighbors are cold channels, the potential for similar fluid heat transfer means the temperature difference is evaluated the same as the hot channels and is summed into a single term. The m in the denominator is the number of discretized steps of the heat exchanger core length. The expanded channels have a rough hexagon shape, which means there are four heat flux terms from neighboring hot channels, and two heat flux terms from neighboring cold channels. Similar fluid heat transfer between cold channels could be in either direction depending on their temperature, and is generally lower in value than the heat transfer with neighboring hot channels. Axial conduction is included in the advanced model with the parameter λ ax , which is a ratio of axial conduction to fluid convection heat transfer.
A w is the cross-sectional area of the HX wall, k w is the thermal conductivity of the HX wall, A f is the cross-sectional area of the HX fluid channel, k f is the thermal conductivity of the fluid, and L HX is the length of the HX core. A factor of half of this is included in the axial conduction term because only one half of the heat conducted in the axial direction enters a channel, since it is conducted through HX walls which have a channel on each side.
Heat leak is included in the advanced model with the parameter λ leak , which is the ratio of external heat leak to the heat transferred within the HX.
The subscript "ins" stands for insulation, which describes the surface area (A), thermal conductivity (k), and thickness (L) of the insulation covering the HX. To derive the factor of 1 4 , consider if the cold inlet temperature of the HX is at environmental temperature. Then the difference between the average HX temperature and the environment is one half of the difference driving heat transfer between fluids, resulting in the factor of 1 4 . Using the overall heat balance: where the flux of each convective side is: . q conv,c,o,j = V C A cs ρ C C p,C T C,j .
The average velocity of fluid is V, the cross-sectional area is A, the fluid density is ρ, the fluid specific heat is C p , and T is fluid temperature.
These equations are solved using MATLAB (The MathWorks, Inc, Natick, MA, USA) for multiple cases, which validate the model. A case with only two channels (no maldistribution or edge effect) compared the advanced model to a simple NTU model, shown in Equation (12): The advanced model exactly reached the NTU result when the number of cells in the model approached infinity. A case with edge effect but no maldistribution was compared to results from a plate heat exchanger by Shah and Focke [53]. Results extrapolated from the advanced model for comparison with a plate heat exchanger were within 0.1% of Shah and Focke. Last, cases with maldistribution and axial conduction compared the advanced model to analytical results using the Peclet number (Pe) and the irregularity parameter (ε): where the Peclet number is the ratio of convective heat transfer to transverse heat transfer as previously defined, and ±a is the channel diameter variation. The advanced model results shared good agreement with this analytical method. Tables 1 and 2 show the default input and output of the MATLAB calculations. See [32] for more details on validation.

Heat Exchanger Prototype
The prototype HX was fabricated out of low-density polyethylene (LDPE) sheets that were 28 µm thick [17] on an open-source polymer laser welding system [15]. The sheets were laser welded together in a pattern that would create both the core and the manifolds of the HX. The HX was expanded with air and fixed into shape. The channels were nominally smooth hexagons with a hydraulic diameter of 2.1 mm in a pattern similar to a honeycomb, shown in Figure 3. The HX was tested with water-to-water flow at a variety of flow rates driven by a water head. The expansion process was not perfect, so a CT scan was used to ascertain the core channel heights as a function of their length [19]. The present paper applies the advanced model to the prototype and discusses the results against experimental data. The conclusions note that this model can be applied to other HXs.

Heat Exchanger Manifold
The manifold was not included in the discretized model because of a lack of CT scan data. The monotonically varying flow rates across the HX in the flow correlation studies [32] could be thought to be due to manifold maldistribution, although the manifold is not explicitly modeled. This section reviews other possible impacts of the manifold for an idealized HX to extract the concepts with the goal of improving the advanced model's fit to the experimental data.
One impact of the manifold is the heat transfer that occurs in it. For the prototype, the heat transfer area in the manifold is approximately 70% the size of the heat transfer area in the core. Therefore, the boundary condition at the core of the cold-water inlet could be significantly hotter, and the hot water inlet could be significantly colder due to heat transfer in the manifold.
For the ideal case of no flow maldistribution, the temperatures of the inlet boundary conditions to the core were varied linearly between two values. Since the cold is entering in the side tube, the first channel will have nearly the same temperature as the tube (see Figure 4). However, the water that goes to the far channel will warm up as it goes through the manifold. The hot inlet is not so straightforward. The path length for the hot flow in the manifold is the same for each channel. However, the side of the HX where the cold flow entered directly into the core will be cooler than the other side of the HX. Therefore, on the cooler side of the HX, the hot flow will cool down more. The extreme case would be the hot flow cooling down as much as the maximum that the cold flow warmed up (symmetric temperature boundary conditions) (see Figure 4). The other extreme case is the hot flows all warming up the same (asymmetric temperature boundary conditions) (see Figure 5). These boundary conditions were used in additional effectiveness trials on an idealized prototype with no channel height variation.
When flow maldistribution is added due to real channel height variations, the situation becomes more complex. The average deviation in temperature from the two boundary inlet conditions is no longer guaranteed to be the same. This is a small problem for randomly generated flow deviations, which were used in previous work for theoretical manifold investigation [32]. However, for the actual prototype channel height, there is some monotonic flow maldistribution. For instance, in order to conserve energy (with equal heat transfer in each manifold), with hot boundary conditions of 0.8 and 0.9, cold boundary conditions of 0 and 0.3 would be expected, but 0 and 0.46 was required. This is not physically realistic, so the solution is to sum up the heat flux of pre-heating the cold water in the manifold and applying that same flux to the hot for post-cooling. Similarly, the heat flux of pre-cooling the hot flow in the other manifold is applied to the cold flow for post-heating.

Flow Maldistribution with Measured Height Data
The measured height of the core of the prototype is shown in Figure 6, split into discretized cells (CT scan data were not available for the manifold). There are 26 cells for each channel in the prototype, and there are 70 channels. The data are organized such that the first 26 cells represent the channel in Cold Row 1-Cold Column 1 (see also Figure 3), the next 26 represent Cold Row 1-Cold Column 2, to the end of Cold Row 1. Then the next channel is Hot Row 1-Hot Column 1, to the end of Hot Row 1 and so on. Starting from Cold Row 1-Column 1, the first cold row ends at cell number 182, which can be seen on Figure 6 as the first sharp decrease in cell height because of the pinched sides on the layer. The next layer starting with Hot Row 1-Hot Column 1 goes from cell 183 to 365, and each layer following in increments of 182 cells. Note the lowest cell height is 0.115 mm which is the pinch height (minimum height required for modeling, see [32]). From these data, the change is more rapid at the top and bottom of the HX (left and right of the figure). The flow predicted by the advanced model in the actual HX based on channel heights in the core is shown in Figure 7, which shows the positive correlation between the hot and cold flows. The irregularity parameter, ε, characterizes the channel flow maldistribution with 1 being no maldistribution [52]. At high NTU, the effectiveness approaches the irregularity parameter. Based on the channel height geometry, the cold ε = 0.591, and the hot ε = 0.540. The sixth and seventh column channels are nearly all pinched off for both cold and hot flow. There is also less flow in the first column than in the center columns. This makes sense because gluing the layers on the sides of the HX pinched the side columns together. The effectiveness results at high NTU are very sensitive to the assumed insulation value for the supply tubes. A combination of fiberglass wrap and foam that was designed to enclose tubes was used, as was feasible for the shape at different points. Results are shown for the reasonable estimate of 1 cm of foam equivalent on average over the tubes.
The 10 experimental trials were simulated in the advanced model. The averaged results of the advanced model are compared with the experimental data in Figure 8. The overall root mean square (RMS) effectiveness error of 3.8% is larger than expected with the effectiveness expected experimental error of 0.7% and the flow rate experimental error of 1-2% (small effectiveness impact). The experimental case had ε = 0.54 and 0.59, but this irregularity parameter increased as high NTU allowed thermal coupling between the channels. A reasonable fit to the data was achieved with rigorous handling of the heat leak. This means that the internal effectiveness of the HX was high (>80%), but the heat leak brought effectiveness down.

Manifold Heat Transfer with Idealized HX
The error between the advanced model and experimental results suggests an unaccounted influence from the manifold. These next sections introduce the results of effectiveness trials at varying inlet temperatures that mimic specific conditions in the manifold, first without channel height variation in the HX.
The cases of NTU = 5 and NTU = 100 correspond to an advanced model effectiveness of 79% and 98%. Since the manifold takes up about 40% of the total heat transfer area, but part of it produces cross flow, which is not as effective, each manifold effectiveness was assumed to be 15%. Therefore, the cold effectiveness in the entry manifold would be 0% and 30% for all cases. The hot effectiveness in the entry manifold would be 0% and 30% for the symmetric case. For the asymmetric case, the hot effectiveness would be 15% and 15%.
The case of NTU = 5 was run first in the advanced model. The case of cold effectiveness in the entry manifold at 15% and 15% and hot effectiveness in the entry manifold at 15% and 15% is not physically realistic for the expanded HX ("equal case"), but it provides an interesting starting point. The overall η = 85.3%, which is consistent with η = 79% from the advanced model between the modified boundary conditions (not shown in graphs). Figures 9-12 show the effectiveness of each channel. The points are connected because the effectivenesses of nearby channels are correlated due to thermal coupling, similar boundary conditions, and similar expansion. Figures 9 and 10 have identical channels (no flow maldistribution). The cold effectivenesses for the symmetric case are shown in Figure 10 (the hot is a mirror image). The cold effectiveness can approach the hot end boundary condition of 70% on the left side and 100% on the right. These graphs correspond to η = 0.679, with an extrapolated value of 0.688. By conservation of energy, if the cold flow was preheated in the manifold 15%, then since these are balanced flow cases, the hot flow would have to be cooled by 15% in the manifold. Therefore, the overall extrapolated η = 0.838. For the asymmetric case the overall extrapolated η = 0.845.    Next the case of NTU = 100 was run. The symmetric case was again mirror images, see Figure 13. The temperature variation between the rows was much less than the NTU = 5 case because there is much stronger thermal coupling in the NTU = 100 case.
As in the NTU = 5 case, the cold effectivenesses for uniform and symmetric flow were approaching the hot boundary conditions of 70% and 100%. The overall extrapolated η = 98.4%.
The NTU = 100 asymmetric case had temperatures that were again not mirror images, and this time the difference was extreme (see Figure 14).  The hot side behaved similarly to both sides of the symmetric case above. There was significant heat transfer from the left to the right side of the HX. The overall extrapolated effectiveness for the asymmetric case is 98.5%, which again is slightly higher than the symmetric effectiveness, again due to lower entropy production.

Manifold Heat Transfer with Measured Geometry
With channel variations from the HX measured geometry, the results from trials with NTU = 100 and symmetric temperature boundary conditions are shown in Figure 15. The extrapolated cold η = 76.8% the extrapolated hot η = 86.0%. However, different amounts of heat transfer occurred in the two manifolds, so the extrapolated overall η = 97.4% for both hot and cold. This is higher than the value of 96.2% without manifold heat transfer, so the manifold heat transfer was still beneficial despite flow maldistribution. For the prototype with NTU = 100 and asymmetric temperature boundary conditions, the cold and hot effectivenesses are shown in Figure 16. The extrapolated cold η = 82.2% and the extrapolated hot η = 85.8%. Again, different amounts of heat transfer occurred in the two manifolds, so the extrapolated overall η = 97.2% for both hot and cold. This is lower than the value of 97.4% for correlated, which is surprising. Perhaps this is due to this specific maldistribution producing more entropy for the uncorrelated case.
For the prototype with NTU = 5 and symmetric temperature boundary conditions, the cold and hot effectivenesses are shown in Figure 11. The extrapolated cold η = 56.2% and the extrapolated hot η = 65.5%. However, different amounts of heat transfer occurred in the two manifolds, so the extrapolated overall η = 76.8% for both hot and cold. This is lower than the value of 79.0% without manifold heat transfer, so surprisingly the manifold heat transfer was detrimental with maldistribution.
For the prototype with NTU = 5 and asymmetric temperature boundary conditions, the cold and hot effectivenesses are shown in Figure 12. The extrapolated cold η = 59.9% and the extrapolated hot η = 63.6%. Again, different amounts of heat transfer occurred in the two manifolds, so the extrapolated overall η = 75.2% for both hot and cold. As for NTU = 100, this is lower than the value of 76.8% for correlated. Manifold heat transfer decreased the overall effectiveness even more in this case.
The conclusion for manifold heat transfer coupled with flow maldistribution for the prototype is that the manifold heat transfer increases effectiveness at high NTU and decreases effectiveness at low NTU.

Discussion
Examining the trials with uniform flow and NTU = 5, Figures 9 and 10, the effectiveness is significantly higher than the similar case with no manifold heat transfer, but lower than the "equal case." This is reasonable because heat transfer between the columns would produce entropy, and the outlet temperatures not being equal cause mixing, producing more entropy. For the asymmetric case, the temperatures are not mirror images, as expected ( Figure 10). The effectiveness is also greater than the symmetric case, which makes sense because there is less entropy production. It is reasonable that the highest effectiveness achieved is the hot left because it is approaching the most extreme temperature of the cold inlet. The cold effectivenesses on both the left and the right are approaching 85% (hot inlet), but the right approaches are closer because it is already preheated in the manifold.
As for the symmetric NTU = 100 cases (Figures 13 and 15), one would expect an approach within 1% of the boundary conditions because of the very high NTU. However, only 93% was achieved on the hot side, but 73% was achieved on the cold side, which is higher than the boundary condition. This means that heat was actually flowing from the "cold" channels to the "hot" channels. That is to say, the flow did not violate the laws of thermodynamics but due to maldistribution the "cold" channels actually became warmer than the "hot" channels. This also explains the very small temperature deviation between rows where this occurred. This is because there is a very small temperature difference between the hot and cold channels, so it did not matter as much that row 1 for cold flow and row 5 for hot flow only have half as many neighbors, because less heat is being transferred regardless of the number of neighbors.
On the cold side of the asymmetric NTU = 100 cases (Figures 14 and 16), one would expect a very close approach (within 1%) to the boundary condition with high NTU, thus yielding near η = 85%. This was true on the right side, but not for the left. The left started out at a colder temperature, but still one could think that the approach should have been closer. This can be explained by the transfer of heat from the left side to the right side of the HX.
The overall conclusion for varying the temperature boundary conditions with no flow maldistribution is that manifold heat transfer increases overall effectiveness roughly what would be expected by the greater heat transfer area in the manifolds.
Another impact of the manifold in addition to altering boundary conditions is the much higher flow resistance for the bent channels. The bent channels are those that enter one side of the HX, turn 90 • to flow axially, and then exit the other side of the HX in a Z type pattern. The straight channels enter one end of the HX and exit the other end of the HX without turning. The relative impact of some of the core channels being pinched would be smaller in this case of bent flow. Therefore, the flow maldistribution on the bent side would be decreased. The current maldistribution on the hot and cold sides is correlated because of the significant pinching on the edges of the HX for both hot and cold. Reducing the flow variation on the bent side would still maintain positive correlation, but it would be less positively correlated. This would reduce the effectiveness somewhat.
A further impact of the manifold is flow maldistribution to the channels. With the Z-type flow, the flow resistance of each channel should be similar. However, there could be significant variations in the flow resistance for different channels in the manifold. This could occur in a uniform way, for instance by having significantly higher flow on one side of the HX than the other. In this case, the heat would have to conduct a much greater distance to smooth out the temperature differences. The Peclet number Pe is the ratio of convective heat transfer to diffusive heat transfer; in this case the latter is the transverse heat transfer. Therefore, the equivalent Pe would be significantly greater with this increased conduction distance.
If heat is conducted from the neighbor of a channel to the neighbor on the opposite side, the Nusselt number should be one (laminar flow). This is about one fourth the case of heating or cooling from all sides of one channel. Therefore, the thermal coupling between channels on opposite sides of the HX would be significantly less. This further reduces the gradient in effectiveness from NTU = 5 to NTU = 100.
It could be possible that the amount of maldistribution changed with flow rate. If maldistribution were significantly higher at high NTU, then the effectiveness could saturate more. A reason for this is the changing pressure experienced by the HX at different flow rates. At 0.01 atm., most of the 100 mm water driving pressure drop at high flow occurred in the manifold for the bent flow. It takes approximately 0.02 atm to inflate the channels to 1 mm high. With some assumptions about the width and length of the pinches, the average height for the bent flow in the manifold was estimated to be approximately 0.2 mm. This is consistent with the observation that the manifold did not contract in a measurable way in the axial direction. Since the pinches in the manifold are only~0.2 mm, a 0.01 atm change in pressure due to flow rate could have significant effects on the maldistribution.

Conclusions
A finite difference model of a HX considered the edge effect, maldistribution, axial conduction, and heat leak. This advanced model is discretized so that it modeled the effects of maldistribution from the first principles, and it has been validated for multiple cases. Data for the channel height as a function of length in the core of the HX were obtained via CT scan. The results run on this core geometry were compared to the experimental results on an expanded microchannel HX. The RMS error was 3.8%.
Though the manifold geometries were not put into the model (because the data were not available), manifold heat transfer was approximated by varying the temperature boundary conditions on the core for the advanced model. With no flow maldistribution, manifold heat transfer increases the overall effectiveness roughly what would be expected by the greater heat transfer area in the manifolds.
However, manifold heat transfer coupled with flow maldistribution for the prototype causes a decrease in the effectiveness at low NTU, and an increase in effectiveness at high NTU.
Another impact of the manifold is the much higher flow resistance for the bent channels. This would reduce the effectiveness somewhat. A further impact of the manifold is flow maldistribution to the channels. The equivalent Pe would be significantly greater, so the effectiveness would saturate more at high NTU.
If heat is conducting from the neighbor of a channel to the neighbor on the opposite side, the Nusselt number should be one if the channels were square. This is approximately one fourth the case of heating or cooling from all sides of one channel. Thus, the thermal coupling between channels on opposite sides of the HX would be significantly less. This further reduces the increase in effectiveness from NTU = 5 to NTU = 100.
It could be possible that the amount of maldistribution changed with flow rate. If maldistribution were significantly higher at high NTU, then the effectiveness could saturate more. Therefore, with further refinements of the advanced model, a better fit to the experimental data may be achieved. If manifold geometry could be obtained, improvements to microchannel heat exchanger modelling could be made by computing manifold heat transfer.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.