Plana Fattori
email: artemio.planafattori@agroparistech.fr

Two-way coupling of fluid-flow, heat-transfer and product transformation during the heat treatment of a starch suspension inside a tubular exchanger

Keywords:

When a liquid product is submitted to thermal treatment under continuous flow, its transport properties can evolve as a consequence of changes occurring on the product constitution (continuous phase, particles...). This two-way coupling is implicitly considered when the residence time distribution (RTD) associated with the product is studied through experimental work. For a number of liquid food products, there is experimental evidence that the minimum residence time decreases as the transformation state progresses. In this study, modeling work is employed for assessing the influence of transformation on the RTD associated with a starch suspension under continuous heat treatment. After including the impact of starch granules swelling on the suspension viscosity, the minimum residence time decreases from 55 % to 41 % of the mean value.

The role played by the transformation in driving the RTD cannot be neglected if the transport properties evolve with the product transformation state.

INTRODUCTION

Many liquid food products are submitted to continuous thermal processing as a way for improving their safety and extending their shelf life [START_REF] Ramaswamy | Residence time distribution (RTD) in aseptic processing of particulate foods: A review[END_REF].

Well-known examples are milk and dairy products, fruit juices, soups, sauces, and egg products. Continuous thermal processing is also employed for inducing transformation of liquid food products, like starch gelatinization [START_REF] Liu | Thermal processing of starch-based polymers[END_REF] and denaturationaggregation of proteins [START_REF] Janhoj | Sensory and rheological characterization of low-fat stirred yoghurt[END_REF][START_REF] Navarra | Thermal aggregation of beta-lactoglobulin in presence of metal ions[END_REF]. Heating (and cooling) sections are employed for rapidly increase (and decrease) the temperature in the liquid stream. Fluid flow, heat transfer and product transformation phenomena are coupled in such cases; indeed, the rheological behavior of a liquid product can evolve as the latter is progressively transformed. These phenomena affect not only the final product characteristics but also the behavior of the velocity field inside the processing unit. This is why the distribution of residence times (associated with the velocity field) can be correlated with the product transformation. The analysis of residence times has been proposed decades ago for studying the performance of chemical processing units [START_REF] Danckwerts | Continuous flow systems: Distribution of residence time[END_REF][START_REF] Levenspiel | Patterns of Flow in Chemical Process Vessels[END_REF]. Since then, the influence of process operating conditions on the final characteristics of liquid food products has been assessed after much experimental work (which is expensive in products, energy and manpower) and modeling efforts (which are hard to be generalized). The experimental techniques used to gather residence time distributions and the mathematical models used to describe them for liquid food products have been reviewed (e.g. [START_REF] Ramaswamy | Residence time distribution (RTD) in aseptic processing of particulate foods: A review[END_REF][START_REF] Torres | Residence time distribution studies in continuous thermal processing of liquid foods: A review[END_REF]. Recent contributions include, for instance, the experimental assessment of the influence of whey protein aggregation on the residence time distribution as measured after heating and holding sections under different flow rates [START_REF] Ndoye | Influence of whey protein aggregation on the residence time distribution in a tubular heat exchanger and a helical holding tube during heat treatment process[END_REF], and the application of computational fluid dynamics (CFD) in modeling fluid flow and heat transfer along a double-pipe helical heat exchanger as a way for better understanding how curvature affects the processing uniformity of product sterilization [START_REF] Rennie | Numerical analysis of the lethality and processing uniformity in a double-pipe helical heat exchanger[END_REF].

The distribution of residence times depends on the apparent viscosity characterizing the liquid product. Since the 1960s [START_REF] Scalzo | Residence times of egg products in holding tubes of egg pasteurizers[END_REF], for egg products) and perhaps earlier, experimental work has put in evidence that an increase of the product viscosity allows a decrease of the minimum residence time. Two limit situations can be identified. The first is associated with relatively low viscosity and isothermal conditions (like milk at 50 °C: about 2 mPa . s); under fully-developed parabolic flow inside a tube, the minimum residence time is then 50 % of the mean residence time. A second limit situation is associated with viscosity high enough for enabling plug flow behavior, like a fully-gelatinized starch suspension: (about 6 Pa . s; [START_REF] Liao | Role of thermo-rheological behaviour in simulation of continuous sterilization of a starch dispersion[END_REF]; the minimum and the mean residence times are then similar because the flat velocity profile. Such limit situations are associated with uniform viscosity fields of liquid products whose characteristics (sizes, etc) are the same in the whole processing unit. Under general cases, the liquid product experience progressive transformation while flowing in the processing unit; it is submitted to heat fluxes across the walls, in whose vicinity fluid parcels often move slower than at the exchanger's axis.

Differences in the velocity can be translated into differences in exposure time to heating or cooling, and hence into different states of product transformation. Therefore we can expect that liquid products characteristics, notably their viscosity, reach non-uniform distributions in the processing unit as a consequence of heat treatment, turning more difficult to predict velocities and residence times.

Let us examine in detail the thermal denaturation-aggregation of whey proteins, which has been studied for a number of reasons including a) the fouling of milk pasteurization units and b) the growing development of dairy products with specific attributes. From the experimental point-of-view, [START_REF] Ndoye | Influence of whey protein aggregation on the residence time distribution in a tubular heat exchanger and a helical holding tube during heat treatment process[END_REF] observed that the occurrence of whey protein aggregation allowed a reduction of the time spent by the fastest fluid parcels. Their Table 2 shows that the minimum residence time in the heating section decreased by about 30 % (from 13.1 to 10.2 s) and in the holding tube by about 10 % (from 44.6 to 40.3 s) after heating the liquid product up to 60 °C (no aggregation) and then to 87 °C (with aggregation) at 49 L/h. Those authors interpreted the reduction of the minimum residence time in terms of an increase in the maximum axial velocity. Such an increase was attributed to the occurrence of aggregation phenomena, which is expected to be more important near the wall because the highest shear rate. In other words, fluid parcels moving near the wall are slowed down while those running at the axis are accelerated, justifying the higher maximum axial velocity and the lower minimum residence time values. Offering a complementary point-of-view from modeling work, [START_REF] Chantoiseau | Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an Eulerian/Lagrangian approach[END_REF] had developed a CFD model for studying the thermal denaturation-aggregation of whey proteins, by coupling fluid flow, heat transfer and particle aggregation. Those authors put in evidence the importance of radial distributions of state variables (temperature, velocity field, particle sizes) in explaining the size distribution of protein aggregates at the outlet. Such a diagnostic tool of the product transformation at every point of the domain is clearly useful in analyzing the process operating conditions and the final characteristics of the liquid product.

The application of mechanistic models constitutes a physically-based strategy for assessing the role played by process operating conditions on residence time distributions. In particular, we advocate the implementation of models coupling fluid flow, heat transfer and product transformation in studying the evolution of liquid food products under heat treatment and continuous flow. Previous efforts have been identified, after neglecting heat transfer while taking into account non-Newtonian twophase behavior [START_REF] Sandeep | Modelling non-Newtonian two-phase flow in conventional and helical-holding tubes[END_REF], or assuming no changes in liquid product viscosity [START_REF] Rennie | Numerical analysis of the lethality and processing uniformity in a double-pipe helical heat exchanger[END_REF][START_REF] Kechichian | Modeling of continuous thermal processing of a non-Newtonian liquid food under diffusive laminar flow in a tubular system[END_REF].

The present study analyzes the influence due to the transformation of a liquid product on the residence time distribution experienced by the own product while it flows along a heat exchanger. Such an influence is estimated with the help of a CFD model coupling the phenomena of fluid flow, heat transfer, and product transformation.

Because the residence time distribution is a consequence of the flow and thermal operating conditions under which the liquid product was transformed, it can be estimated as a post-processing task after running the coupled model. In other words, the residence time distribution can be estimated from the velocity field resulting from the numerical simulation. Once the CFD model couples the phenomena of interest, the velocity field incorporates the effects due to the transformation on the fluid viscosity.

In order to put in evidence the two-way coupling between fluid flow, heat transfer and product transformation as well as the importance of such a coupling on the residence time distribution, we retained a liquid food product whose progressive evolution can allow significant consequences on the fluid viscosity: an aqueous suspension of starch granules. Under heat treatment in excess water, starch granules can experience phase transitions whose extent can depend on the heating rate and on the final temperature reached. The concept of "gelatinization" refers to the destruction of the crystalline structure in granules; such an irreversible process includes, in a broad sense and in time sequence, the granule swelling, the melting of native crystalline structure and the molecular solubilization [START_REF] Liu | Thermal processing of starch-based polymers[END_REF]. We focus our attention on the granule swelling, whose consequences on the viscosity of starch suspensions have been observed at least since the 1940s [START_REF] French | Physical properties of starch[END_REF]; within a small range of temperature at about 65 °C, starch granules suddenly increase many times in size while the viscosity rapidly rises. More recently, it has been shown that the viscosity of the starch suspension relative to that of water depends not only on the volume occupied by the granules but also on their shape and deformability [START_REF] Ellis | A comparison of the viscous behaviour of wheat and maize starch pastes[END_REF]. Starch research is a very active subject, and much effort has been done in characterizing native and modified starch, including various genotypes and climates (see for instance the review by [START_REF] Singh | Morphological, thermal and rheological properties of starches from different botanical sources[END_REF].

The evolution of a starch suspension in hereafter studied with the help of numerical simulations, whose results are later analyzed in terms of the residence time distribution. Next section describes the methods employed for solving the coupled problem of fluid flow, heat transfer, and starch granule swelling. The case study involves a starch suspension for which a thermo-rheological model has been built through experimental work. The following section presents and discusses the more relevant results regarding: a) the evolution of state variables (temperature, velocity field, and swelling degree) along the domain, b) the importance of thermal histories on the transformation state reached by the suspension, and c) the residence time distribution estimated from the velocity field obtained under selected scenarios. The final section summarizes the work accomplished and emphasizes our conclusions.

METHODS

The evolution of a liquid product under heat treatment and continuous flow is hereafter studied with the help of numerical simulations. For an incompressible fluid, the conservation equations for mass, momentum and energy under steady-state conditions can be written as [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF]:

 . u 0 (1)  u .u  . ( p I  u (u) T ) ) (2)   CP u .T  . ( kT ) (3)
where state variables are the velocity vector u (m/s), the pressure p (Pa), and the temperature T (K); bulk properties of the liquid product are the density  (kg/m 3 ), the apparent viscosity  (Pa . s), the specific heat capacity at constant pressure CP (J/(kg . K)), and the thermal conductivity k (W/(m . K)).

Equation ( 2) requires the apparent viscosity  associated with the liquid product under consideration, in this case an aqueous suspension of starch granules. For a suspension of rigid particles, the suspension viscosity  can be estimated by multiplying the continuous phase (pure water) viscosity W and the relative viscosity R [START_REF] Bird | Transport Phenomena, 2nd Edition[END_REF]. The increase of relative viscosity with the volume fraction  occupied by the solid particles in the suspension has been studied for a number of systems (see for instance Chen and Law (2003) for a review). [START_REF] Barnes | An Introduction to Rheology[END_REF] have proposed two correlations for the increase of R with  under two limit situations, associated with very low and high shear rate respectively. Further, those authors suggested that suspensions associated with  < 0.3 could be assumed Newtonian, because those two limit correlations allow essentially the same viscosity. Under higher solid volume fractions, the rheological behavior of starch suspensions may not only differ from the Newtonian one but also depend on parameters like the granule deformability and the starch solubility [START_REF] Doublier | A rheological investigation of cereal starch pastes and gels. Effect of pasting procedures[END_REF]; see also [START_REF] Lagarrigue | The rheology of starch dispersions at high temperatures and high shear rates: A review[END_REF] for a review). In summary, reliable estimates of relative viscosity can be obtained after some knowledge of the own state of transformation of the starch suspension.

The transformation state of a starch suspension can be characterized by the swelling degree of starch granules in water [START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF]:

S = (D -D0) / (DM -D0), ( 4
)
where D is the mean diameter of starch granules, and D0 and DM are its values before heat treatment (native starch) and after complete treatment, respectively. Mean volume-weighted diameters are here employed because the transformation of the liquid product is associated with progressive swelling of starch granules and consequent increase of their volume fraction in the aqueous suspension. Such an increase can be translated into relative viscosity (see below). The swelling degree is a continuous function assuming values from 0 (native granules only) and 1 (maximum swelling condition). Its variation with time can reasonably be represented through a second-order kinetics equation [START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF]:

dS / dt = V ( 1  S ) 2 , ( 5
)
where V is the reaction rate constant which is estimated through the Arrhenius law,

V(T)V0 exp(EA / (RT) ) . (6)
Under steady-state conditions and after including a diffusion term for starch granules in the suspension, equation ( 5) can be written as a convection-diffusion equation:

u .S  V ( 1  S ) 2 +. (D SS ), (7)
where D S is the diffusion coefficient (m 2 /s).

In the present study, moderate heating conditions are applied to a given starch suspension. As indicated by experimental work; the maximum volume fraction of the swollen granules is about 0.5 (see below). Under such conditions, the relative viscosity depends mostly on the solid volume fraction [START_REF] Lagarrigue | The rheology of starch dispersions at high temperatures and high shear rates: A review[END_REF]. The model proposed by [START_REF] Thomas | Transport characteristics of suspension: VIII. A note on the viscosity of Newtonian suspensions of uniform spherical particles[END_REF],

R = 1.510.05  2  0.00273 exp( . (8)
is hereafter employed, because it provides an increase of R with  which falls within the interval bounded by the two limit correlations proposed by [START_REF] Barnes | An Introduction to Rheology[END_REF] (see for instance the Figure 3 of Chen and Law ( 2003)). The solid volume fraction can be evaluated from the mean volume-weighted diameter:

0DD0 3 , (9)
where 0 and D0 values characterize the starch suspension before any treatment, while D corresponds to the present local state of the starch suspension:

D = D0 + S ( DM -D0 ) . ( 10
)
Much effort has been accomplished in characterizing native and modified starch, including different genotypes and various climates [START_REF] Yeh | Kinetics of phase transition of native, cross-linked, and hydroxypropylated rice starches[END_REF][START_REF] Buleon | Starch granules: Structure and biosynthesis[END_REF][START_REF] Paterson | Rheology and granule size distributions of corn starch dispersions from two genotypes and grown in four regions[END_REF][START_REF] Nayouf | Effect of thermomechanical treatment on the rheological properties of crosslinked waxy corn starch[END_REF][START_REF] Singh | Morphological, thermal and rheological properties of starches from different botanical sources[END_REF][START_REF] Choi | Swelling characteristics of native and chemically modified wheat starches as a function of heating temperature and time[END_REF][START_REF] Lindeboom | Analytical, biochemical and physicochemical aspects of starch granule size, with emphasis on small granule starches: A review[END_REF][START_REF] Xie | Rheological properties of starches with different amylose/amylopectin ratios[END_REF][START_REF] Hoover | The impact of heat-moisture treatment on molecular structures and properties of starches isolated from different botanical sources[END_REF]. Equations (1, 2, 3, 7) were solved under steady-state conditions with the finiteelement method (simulation package COMSOL Multi-physics, version 3.5a; see [START_REF] Zimmerman | Multiphysics Modelling with Finite Element Methods[END_REF]). Second-order Lagrange finite elements were employed in all cases, except for the pressure for which first-order finite elements were assumed. We choose the Parallel Sparse Direct Linear Solver (PARDISO) for solving the large systems resulting from the linearization of the coupled equations under consideration.

PARDISO implements an efficient method for large systems of equations like those arising from the finite element method [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. Results discussed below were obtained by assuming a relative tolerance of 10 6 . No apparent reasons forbid the solution of this coupled problem with the help of other modeling software or programming language.

In our case study, the gradients of velocity, temperature, swelling degree and suspension viscosity are mainly oriented along the domain's radius; therefore dr has to be small. The numerical representation of the coupled processes of interest is difficult at the heating wall because of the relatively high magnitude of velocity and temperature gradients that prevail there; both swelling degree and suspension viscosity are expected to increase sharply near the corner of the domain inlet with the heating wall. The representation of such a challenging behavior can be achieved with the help of equilateral triangles or squares. The geometry of the computational domain suggests the adoption of meshes constituted of rectangular cells, characterized by dimensions dz and dr along the domain's length (Z = 1 m) and radius (R = 5 mm), respectively. To mesh the domain into a great number of squares could be a suitable strategy for representing the case study of interest. However, calculator memory capabilities can limit the overall number of squares to be considered.

Sensitivity tests were conducted in order to assess the influence of the mesh strategy (squares or rectangles) and mesh resolution (as indicated by the number of cells) on selected results. Figure 2 presents mass-weighted (bulk) estimates of temperature (left) and starch swelling degree (right) at the domain outlet as functions of the number of cells; results were obtained after assuming selected values of the aspect ratio dz / dr: 1 (squares), 5 ("short rectangles") and 10 ("long rectangles"). Bulk estimates exhibit a near asymptotic behavior as the number of cells increases; such a behavior appears when partial differential equations are solved through numerical methods with the help of grids progressively finer. Beyond 200,000 cells, the application of short rectangles provides essentially the same bulk results than employing squares, while the application of long rectangles provides slightly lower values. The calculation time for the converged solution depends on the mesh resolution, from few minutes for 50,000 cells up to about 20 hours for 250,000 cells. It must be noted that, for a given number of cells, gradients along the domain's radius are better represented by employing long rectangles; further, gradients along the domain's length are better represented by employing squares.

Equation ( 7) requires the specification of the diffusion coefficient D S, which can be assessed through the analysis of the Brownian motion in the case of a dilute suspension of spherical colloidal particles [START_REF] Bird | Transport Phenomena, 2nd Edition[END_REF]. The diffusion coefficient depends on the continuous phase temperature and on the particle size; in our case study D S can vary along the evolution of the liquid product in the heat exchanger. A rough estimate can be reached by assuming native starch granules (see above: DO = 15.5 m)

at the inlet temperature (20 °C): we find D S ~ 3 10 -14 m 2 /s. As indicated by sensitivity tests (not shown here), the application of such a value could allow model convergence only after considering a mesh resolution beyond our present memory capabilities (a 64bits calculator disposing of 24 Gb of RAM).

Looking for the numerical solution of the coupled problem, we adopted the highest mesh resolution with short rectangles (dz / dr = 5) which is allowed by our memory capabilities, namely the mesh constituted of 80 rectangular cells along the domain's radius and 3200 along the domain's length (total number of cells: 80 x 3200 = 256,000). The application of short rectangles represents a compromise, because it provides a finer representation of gradients along the domain's radius than with squares, and a finer representation of gradients along the domain's length than with longer rectangles. At such mesh resolution, model convergence was reached after assuming the diffusion coefficient D S = 10 12 m 2 /s in equation ( 7).

RESULTS AND DISCUSSION

The evolution of an aqueous suspension of starch granules along its progression in the heat exchanger is hereafter discussed from model results. Variations along the domain's radius are expected, as consequence of the roles played by the heating wall on phenomena of fluid flow and of heat transfer, and hence of product transformation.

Variations along the domain's length are also expected, because radial distributions of key variables (temperature, velocity field, swelling degree) evolve while the product progresses in the exchanger. finally, that all starch granules present in the suspension are on native condition (that is, null swelling degree). Top-right display of Figure 3 shows that, at the heating wall, the maximum heat treatment (S = 1) takes place in the first centimeters of the exchanger.

The swelling degree decreases with the distance from the heating wall, and such a decrease tendency becomes weaker as the product progresses in the exchanger. There is a weak increase tendency of the swelling degree at the axis of symmetry, consequence of the progressive warming of this region.

The starch dispersion was assumed to contain native granules only at the domain's inlet (diameter = D0 = 15.5 m), i.e. it was associated with the minimum solid volume fraction and relative viscosity values (0  0.03 and R ~ 1.09); as a result, the suspension viscosity at the domain's inlet approaches that for water (about 1 mPa . s).

Bottom displays of Figure 3 show that, with the progressive swelling of starch granules, their mean diameter D and their volume fraction  increase. The impact of granule swelling in increasing the relative viscosity is quite explosive near the heating wall (bottom-center display): the relative viscosity reaches a factor higher than 13 in the first centimeters of the exchanger. Such an increase explains in a large extent the behavior of the suspension viscosity near the heating wall. The representation of such a behavior is very challenging from the numerical point-of-view because the importance of gradients of the variables of interest along both domain's radius and length.

The suspension viscosity depends not only on the volume fraction occupied by starch granules but also on the temperature, via the water viscosity. As shown in the bottom-right display of Figure 3, the suspension viscosity can even decrease with respect to its value at the domain's inlet at regions where the swelling degree is relatively weak (say, below 0.2). Along these regions, the increase of relative viscosity (due to granule swelling) is more than counterbalanced by the decrease of water viscosity (due to temperature increase). An interesting feature can be appreciated in the bottom-right display of Figure 3: because the role played by the temperature influence on the water viscosity, the maximum suspension viscosity value does not occurs at the domain's outlet. This is because the relative viscosity has already reached its maximum, but water viscosity is still decreasing with increasing temperature.

Near the wall, the overall result is the increase of suspension viscosity, due to the increase of solid volume fraction. Such an increase slows down the velocity near the wall. For instance, at r = 0.004 mm (i.e. 1 mm from the heating wall), from z = 25 to z = 75 cm, the suspension viscosity more than doubles (from ~0.75 to ~1.6 mPa . s) while the velocity magnitude decreases by a factor of about two (from ~0.025 to ~0.012 m/s). The evolution of liquid product properties can also be studied with the help of Lagrangian trajectories associated with selected fluid parcels released at the domain inlet. Each fluid parcel can be imagined as a suspension droplet containing an ensemble of starch granules at similar thermodynamic and kinetic conditions. In its right section, More generally, each fluid parcel of a liquid product takes a different route and spends a different time to flow through the processing unit. The distribution of these times at the outlet characterizes the residence time distribution of the process. Since [START_REF] Danckwerts | Continuous flow systems: Distribution of residence time[END_REF], the analysis of residence time distributions has been usefully applied to study the flow pattern, to determine the degree of mixing and to diagnose flow problems such as recirculation, channeling, short-circuiting or stagnation [START_REF] Gutierrez | Residence time distribution in holding tubes using generalized convection model and numerical convolution for nonideal tracer detection[END_REF].

In this study the residence time distribution is estimated from velocity numerical simulations. The residence time distribution is often assessed at the processing unit's outlet; its estimation is here performed by backward integrating the velocity field, from the domain's outlet to the inlet. Such a procedure was implemented for three scenarios To our knowledge, no experimental studies have compared the residence time distribution with and without product transformation in the case of starch suspensions.

However, we have identified two studies about residence time distributions associated with liquid food products whose viscosity increases as the transformation progresses:

the thermal denaturation-aggregation of whey proteins inside a tubular heat exchanger [START_REF] Ndoye | Influence of whey protein aggregation on the residence time distribution in a tubular heat exchanger and a helical holding tube during heat treatment process[END_REF], and the sorbet freezing inside a surface scrapped heat exchanger [START_REF] Arellano | Experimental study and modelling of the residence time distribution in a scraped surface heat exchanger during sorbet freezing[END_REF]. These studies put in evidence that the minimum residence time decreases when the transformation degree becomes higher, i.e. when the heat treatment is assigned to be stronger in the case of whey proteins denaturation-aggregation, and when the wall temperature is assigned to be lower in the case of sorbet freezing.

SUMMARY

Looking for a better understanding of the influence due to the transformation of a liquid product on the residence time distribution experienced by the own product while it flows along a heat exchanger, we estimated residence time distributions from the velocity field resulting from numerical coupled simulation. Such a post-processing task can be performed not only at the domain's outlet, where RTDs is often accessible through experiment, but also at intermediate positions of the domain, what can be particularly useful in design and optimization of processing units.

The approach was illustrated in studying the evolution of a starch suspension in a heat exchanger. Quite standard conditions were assumed (cylindrical geometry, laminar Newtonian flow, uniform wall heat flux, steady-state conditions), and a kinetics equation was adopted for evaluating the swelling degree of starch granules in the aqueous suspension. Despite such an apparent simplicity, the problem is challenging.

Firstly, velocity and temperature gradients are higher in the vicinity of the heating wall.

Secondly, starch granule swelling depends on thermal histories, which are progressively longer towards the wall. Finally, swelling severely impacts the suspension viscosity in the vicinity of the wall and more particularly near the corner with the domain's inlet. In order to satisfactorily represent the coupled processes under consideration, the combination of these three issues strongly suggested us to build a mesh as fine as possible, within the limits presented by our computational capabilities. Our results can be verified by others through any modeling software able to solve a coupled problem involving fluid flow, heat transfer and a kinetics reaction.

Our results put in evidence the importance of separating a) the influence of temperature on the continuous phase viscosity from b) the influence of transformation on the relative viscosity associated with the solid fraction of the suspension; their effects can be opposite. The influence of temperature on the continuous phase viscosity depends on the present local temperature of the liquid product, and in this sense this effect can be considered instantaneous. On the other hand, the influence due to the transformation on the relative viscosity depends on the entire thermal history of fluid parcels. This second effect can be interpreted as the result of all the previous states experienced by the liquid product under the operating conditions of interest.

We advocate the use of a numerical model able to couple the phenomena of fluid flow, heat transfer, and liquid food product transformation in this kind of problem;

otherwise, the velocity field assumed in RTD estimates can be unrealistic. As shown in Figure 5, the minimum residence time can be overestimated by about 30 % (that is, 55% instead 41 % of the mean residence time) if transformation processes are not represented. Product transformation processes have to be represented in order to enable consistent RTD estimates, as suggested by experimental work in similar fields of food engineering, which proved that RTD can significantly depend on transformation occurring inside heat exchangers.

Residence time distributions can be estimated from velocity fields resulting from numerical coupled simulation of any sort, including the representation of transformation processes as complex as needed. In this study we focused the attention on a simple transformation process (swelling of starch granules in water, under continuous heating) whose representation can be achieved through a simple methodology (a second-order kinetics equation). Such a choice was deliberate: looking for putting in evidence the two-way coupling between fluid flow, heat transfer and product transformation, we choose a liquid product whose progressive evolution allows significant consequences on the fluid viscosity. Work in progress includes the coupling of fluid flow and heat transfer with a more complex transformation process found in liquid food engineering, namely the particle aggregation. Its representation involves the population balance equation, and preliminary results are available [START_REF] Chantoiseau | Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an Eulerian/Lagrangian approach[END_REF].

  [START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF] have studied selected starch suspensions, following the kinetics of granule swelling when applying thermo-mechanical treatments in a laboratory simulator coupled to a rheometer. Data on granule swelling and apparent viscosity development were obtained under controlled temperature and shear conditions. In this study we focused the attention on waxy maize starch; its solubility under heat treatment is negligible because the low level of amylose molecules in its constitution. Further, we restrict the analysis to waxy maize starch which has been chemically modified, preventing the rupture of swollen granules. Hereafter, we consider the evolution of the 3.1 % w/w modified waxy maize starch suspension studied by[START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF], who obtained D0 15.5 and DM  39.6 m, 0  0.03, V0 2.9 10 12 s  , and EA  84.8 kJ/mol. The maximum swelling condition (S = 1) corresponds to the solid volume fraction 0DMD0 3 ~ 0.5, which provides an increase of relative viscosity R of about 16. Such an increase cannot be directly translated into the suspension viscosity , because the water viscosity W decreases with the suspension temperature. The influence of temperature on water viscosity can be assessed with the help of standard tables (e.g.[START_REF] Iapws | Supplementary Release on Properties of Liquid Water at 0.1 MPa[END_REF].Coupled phenomena of fluid flow, heat transfer and starch granule swelling are illustrated inside a tubular heat exchanger with length 1 m and radius 5 mm, represented through a two-dimension axial-symmetric domain. In equations (2) and (3), the density , the specific heat capacity at constant pressure CP, and the thermal conductivity k can be estimated by averaging the values associated with water and starch granules, where the mass fraction plays the role of weighting factor (e.g.Delgado et al., 2006). These properties are hereafter assumed the same as for water because the relatively low mass fraction associated with starch granules in the suspension here considered (3.1 % w/w, see above). A uniform flux of 10 kW/m 2 was applied on the heating wall, and the temperature of the starch dispersion was assumed to be 20 °C at the inlet. A fully developed flow profile was assumed at the inlet, and the flow rate was fixed at 10 liters per hour (mean velocity 0.035 m/s). According to the Reynolds number (about 350 at the domain inlet), the flow regime is laminar. Such boundary conditions were chosen in order to allow temperatures up to about 90 °C, which corresponds to a relevant degree of starch granule swelling according to[START_REF] Lagarrigue | Swelling kinetics of waxy maize and maize starches at high temperatures and heating rates[END_REF]. Geometry and boundary conditions are summarized in Figure1.

Figure 3

 3 Figure 3 presents the radial distribution of selected variables at three distances from the domain's inlet: 25, 50, and 75 cm (that is, at 25, 50, and 75 % of the domain's length). Top displays show state variables (temperature, magnitude of velocity field, and swelling degree) obtained by solving equations (1, 2, 3, 7) as explained above. Bottom displays show variables obtained from the former: a) the mean volumeweighted diameter from the swelling degree by applying equation (10), b) the relative viscosity by applying equation (8) after estimating the solid volume fraction through equation (9); and c) the suspension viscosity by multiplying the relative viscosity and the water viscosity at the temperature under consideration.Inlet boundary conditions include uniform temperature (T = 20 °C). Top-left display of Figure3shows that, along the exchanger, temperature increases faster near the wall (r = 0.005 m) than at the axis of symmetry (r = 0); further, the thermal boundary layer, along which the temperature decreases linearly with the distance from the heating wall, becomes progressively wider. Inlet boundary conditions include also a parabolic profile for the velocity field (fully developed flow). Top-center display of Figure3shows that the shape of such a profile evolves along the exchanger: fluid parcels travelling in the vicinity of the wall are slowed down, while those travelling at

Figure 4

 4 Figure 4 summarizes, in its left section, the whole fields of temperature, mean diameter and suspension viscosity. The suspension viscosity field was deformed according to the velocity profiles at the domain's inlet and outlet. Vertical dotted lines indicate the three axial coordinates corresponding to the radial distributions shown in Figure 3. These results put in evidence the evolution of state variables along the heat exchanger, under heating conditions which can be seen as moderate when compared with those assumed by other authors. For instance, Liao et al. (2000) have modeled the evolution of 4 % w/w waxy rice starch suspensions under different operating conditions (exchanger diameter, heating wall temperature, and fluid volumetric flow rate). As a consequence of their heating conditions (up to 145 °C), the increase of suspension viscosity was important enough for causing velocity profiles almost flat (plug flow behavior) where the suspension was very viscous. This was not observed in the present study. Our results correspond to weaker heating conditions (up to 91 °C near the wall) and meaningfully weaker suspension viscosity values (up to 7 mPa . s, instead 6 Pa . s obtained by Liao et al.); as shown in the top-center display of our Figure 3, the resulting velocity field is slightly influenced by the increase of suspension viscosity, while remaining far from the plug flow behavior.

Figure 4

 4 Figure4displays selected properties associated with fluid parcels released at 1 and at 0.5 mm from the heating wall; they spent about 36 and 86 s respectively before reaching the outlet. Fluid parcels released the nearest to the heating wall are associated with higher values of temperature and suspension viscosity. The suspension viscosity also depends on the continuous phase (water) viscosity, which weakly decreases with

  Figure5, the minimum residence time is equal to 50 % of the mean value. The second

FIGURE CAPTIONS Figure 1 :

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: Bulk temperature (left) and swelling degree (right) values at the domain outlet as functions of the total number of mesh cells, as obtained with the CFD model developed in this study. Three sets of results are shown, corresponding to selectedvalues of the aspect ratio dz / dr: 1 (squares), 5 ("short rectangles") and 10 ("long rectangles"). Dimension dz (dr) is along the domain's length (radius) of every rectangular cell constituting the mesh considered. Dotted lines were included for illustration purposes only.

Figure 3 :

 3 Figure 3: Radial distributions of selected variables at selected distances from the domain's inlet (z = 25, 50, and 75 cm).

Figure 4 :

 4 Figure 4: Selected results provided by modeling the coupled problem. Left section presents the temperature (top), swelling degree (center) and suspension viscosity (bottom) fields in the domain under consideration; inlet is on the left, heating wall is on the bottom. Black and white colors indicate minimum and maximum values: about 20 and 91 °C for temperature, zero and 0.9996 for swelling degree, and 0.69 and 7.1 mPa . s for suspension viscosity, respectively; dotted lines indicate the three axial coordinates which correspond to the radial distributions shown in Figure 3. Right section presentsthe time series of selected properties corresponding to the trajectories followed by two hypothetical fluid parcels released at the domain inlet at 0.5 and 1 mm from the wall.

Figure 5 :

 5 Figure 5: Cumulative residence time distributions resulting from the coupled model including fluid flow, heat transfer and starch swelling (white squares), and from successively simpler versions that discard starch swelling (grey symbols) and discard both heat transfer and starch swelling (black symbols).

ACKNOWLEDGEMENTS

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7 / 2007-2013) under the grant agreement n° FP7-222 654-DREAM.